{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第四章　马尔可夫链蒙特卡洛方法（Markov Chain Monte Carlo, MCMC）\n",
    "\n",
    "\n",
    "作者：[王何宇](http://person.zju.edu.cn/wangheyu)\n",
    "\n",
    "[浙江大学数学科学学院](http://www.math.zju.edu.cn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们现在考虑如何为一个随时间变化的系统建模。首先还是需要一些基本假设。\n",
    "## 马尔可夫链\n",
    "假设一个随时间变化的系统，只有有限个状态，分别记作\n",
    "$$\n",
    "\\Omega = \\{x_1, x_2, \\cdots, x_n\\}.\n",
    "$$\n",
    "然后这里时间也是离散整点的，于是我们需要考虑的就是在 $t = 1, 2, \\cdots$ 时刻，系统的状态 $X_t \\in \\Omega$。如果，$X_t$ 的取值仅由 $X_{t - 1}$ 决定，和更早时刻的系统状态无关（**遗忘性**），则称满足此要求的状态序列 $\\{X_t\\}, t = 0, 1, \\cdots$ 为一个**马尔可夫链**。这里 $X_t$ 故意大写，是为了表示它是一个不确定的随机值，它可能是 $\\Omega$ 中任何一种状态，它实际上表示的是在 $\\Omega$ 中取值的一种随机分布，$\\Omega$ 中每一种状态，都以一个确定（仅在 $t$ 时刻）的概率取到。严格地说，马尔可夫链是一种随机分布序列，$t$ 时刻的分布 $X_t$ 仅由上一时刻 $X_{t - 1}$ 的分布决定。因此，如果我们确定给一个 $t - 1$ 时刻的系统状态，也即 $X_{t - 1}$ 的值，比如\n",
    "$$\n",
    "X_{t - 1} = x_i \\in \\Omega,\n",
    "$$\n",
    "那么，$X_t$ 的**分布**就只能有这个 $x_i$ 状态所**完全**确定。\n",
    "记\n",
    "$$\n",
    "P\\{x_i \\to x_j\\} := p_{ij}, \\quad i, j = 1, 2, \\cdots, n.\n",
    "$$\n",
    "则显然有\n",
    "$$\n",
    "\\sum_{j = 1}^n p_{ij} = 1, \\quad i = 1, 2, \\cdots, n.\n",
    "$$\n",
    "（这就是“完全”的意义。）\n",
    "\n",
    "形式上，这个矩阵\n",
    "$$\n",
    "P := \\left[p_{ij}\\right]_{n \\times n}, \\quad i, j = 1, 2, \\cdots, n.\n",
    "$$\n",
    "包含了一个马尔可夫链的全部信息，它的第 $i$ 行，代表了状态 $x_i$ 在下一个时刻可以转变的全部**可能**状态的分布，这里我们指出：\n",
    "1. 对角元 $p_{ii}$ 可以大于零，也即存在一定的概率，在下一个时间系统保持状态 $x_i$ 不变；\n",
    "2. 允许 $p_{ij} = 0$，也即存在这样的状态 $x_j$，它不可能由 $x_i$ 在一个时间步转变而来；\n",
    "3. $p_{ij} > 0$，$i, j = 1, 2, \\cdots, n$.\n",
    "\n",
    "注意这里我们实际上约定了 $p_{ij}$ 是一个和时间 $t$ 无关的常量。这样的马尔可夫链，称为齐次的（homogenous）。在接下去的讨论中，我们将模型局限在有限、齐次马尔可夫链下。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例** 令 $\\Omega = \\{x_1, x_2, x_3, x_4\\}$，转移矩阵为\n",
    "$$\n",
    "\\mathbf{P} = \\left[\n",
    "\\begin{array}{cccc}\n",
    "p_{11} & p_{12} & 0 & p_{14} \\\\\n",
    "p_{21} & 0 & p_{23} & 0 \\\\\n",
    "p_{31} & 0 & 0 & p_{34} \\\\\n",
    "0 & 0 & p_{43} & p_{44}\n",
    "\\end{array}\n",
    "\\right].\n",
    "$$\n",
    "或者等价地，可以用转移图表示为：\n",
    "<img src=\"graph.pdf\" width=\"25%\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们假定初始状态为\n",
    "$$\n",
    "X_0 = x_1,\n",
    "$$\n",
    "则初始状态是确定的，也即\n",
    "$$\n",
    "\\mathbf{p}_0 = (1, 0, 0, 0).\n",
    "$$\n",
    "而下一个状态 $X_1 = x_j$ 的概率分布是 $p_{1j}$，也即\n",
    "$$\n",
    "P\\{X_1 = x_j\\} = p_{1j}, \\quad j = 1, 2, 3, 4.\n",
    "$$\n",
    "\n",
    "一般地，若用\n",
    "$$\n",
    "\\mathbf{p}_t = (p_1, p_2, p_3, p_4)\n",
    "$$\n",
    "来表示 $t$ 时刻系统状态变量 $X_t$ 的概率质量分布（PMF），即\n",
    "$$\n",
    "\\mathbf{p}_j = P\\{X_t = x_j\\}, \\quad j = 1, 2, 3, 4,\n",
    "$$\n",
    "则称向量 $\\mathbf{p}$ 为概率向量（一般约定，它是个行向量）。由马尔可夫链定义，知\n",
    "$$\n",
    "\\mathbf{p}_{t + 1} = \\mathbf{p}_t \\mathbf{P}.\n",
    "$$\n",
    "注意这里行向量用法，进而\n",
    "$$\n",
    "\\mathbf{p}_{t + 1} = \\mathbf{p}_t\\mathbf{P} = \\mathbf{p}_{t - 1}\\mathbf{P}^2 = \\cdots = \\mathbf{p}_0 \\mathbf{P}^{t - 1}.\n",
    "$$\n",
    "到此我们可以给出结论：一个马尔可夫链的演化过程，完全由系统初始状态和转移矩阵确定。由上一章所学，我们可以在一个实际模拟中，通过随机抽取，来模拟一个马尔可夫链。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例** 假设股市是一个有限齐次马尔可夫链（韭菜田），并记三种状态分别为涨、平、跌对应 $x_1$, $x_2$ 和 $x_3$，其转移矩阵为\n",
    "$$\n",
    "\\mathbf{P} = \\left[\n",
    "\\begin{array}{ccc}\n",
    "0.3 & 0.2 & 0.5 \\\\\n",
    "0.4 & 0.2 & 0.4 \\\\\n",
    "0.4 & 0.3 & 0.3\n",
    "\\end{array}\n",
    "\\right],\n",
    "$$\n",
    "于是只要有一个初始状态，我们就可以对股市进行模拟。比如，记初始状态为涨，即\n",
    "$$\n",
    "\\mathbf{p}_0 = (1, 0, 0),\n",
    "$$\n",
    "则经过充分长时间的变化后，会发生什么？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "# 这句话是非标准的python，用于ipthon或jupyter这样的系统中，表示绘图即刻自动展开。\n",
    "%matplotlib inline\n",
    "\n",
    "# 这里把全部Warning过滤掉了. \n",
    "# 参见https://docs.python.org/2/library/warnings.html\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "#np.random.seed(250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们设置初值和转移矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = np.array([1, 0, 0])    \n",
    "P = np.array([[0.3, 0.2, 0.5], [0.4, 0.2, 0.4], [0.4, 0.3, 0.3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后每一次状态变换，都是一次向量乘以矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3, 0.2, 0.5])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = p @ P\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "重复上述操作多次，最后我们发现，\n",
    "$$\n",
    "\\mathbf{p}_t \\to (0.36363636, 0.23966942, 0.39669421), \\quad t \\to \\infty.\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一般地，如果马尔可夫链链的概率向量 $\\mathbf{p}_t$ 当 $t \\to \\infty$ 时趋于一个常向量，我们称其为马尔可夫链的**不变分布**，或**平稳分布**，或**均衡分布**。自然界中有很多现象，是满足这个模型的，但是股市显然不是。\n",
    "\n",
    "从不变分布的定义，我们可以看到，事实上，不变分布就是迭代过程\n",
    "$$\n",
    "\\mathbf{p}_{t + 1} = \\mathbf{p}_t \\mathbf{P}\n",
    "$$\n",
    "的不动点 $\\pi$，满足\n",
    "$$\n",
    "\\pi \\mathbf{P} = \\pi,\n",
    "$$\n",
    "或者就是转移矩阵 $\\mathbf{P}$ 的特征值为 $1$ 对应的左特征向量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理（Perron-Frobenius）** 若转移矩阵 $\\mathbf{P}$ 满足\n",
    "$$\n",
    "p_{ij} > \\delta > 0, \\quad i, j = 1, 2, \\cdots, n, \n",
    "$$\n",
    "则\n",
    "1. $\\mathbf{P}$ 存在特征值 $1$，且对应的左特征向量 $\\mathbf{w}$ 全部分量严格为正；\n",
    "2. 若 $\\mathbf{w}$ 归一化，则\n",
    "$$\n",
    "\\lim_{t \\to \\infty}\\mathbf{P}^t = \\mathbf{1w};\n",
    "$$\n",
    "3. 对于任何概率（行）向量 $\\alpha$，有\n",
    "$$\n",
    "\\lim_{t \\to \\infty} \\alpha \\mathbf{P}^t = \\mathbf{w},\n",
    "$$\n",
    "这里 $\\mathbf{w}$ 是已经归一的，$\\mathbf{1}$ 是分量均为 $1$ 的列向量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定理证明略，我们这里简单说明一下为何 $\\mathbf{P}$ 必然有特征值 $1$。事实上，令 $\\lambda$ 是 $\\mathbf{P}$ 的特征值，$\\mathbf{w} = (w_1, w_2, \\cdots, w_3)$ 是其左特征向量，即：\n",
    "$$\n",
    "\\mathbf{w} \\mathbf{P} = (\\sum_{i = 1}^nw_ip_{i1}, \\sum_{i = 1}^nw_ip_{i2}, \\cdots, \\sum_{i = 1}^nw_ip_{in}) = \\lambda \\mathbf{w} = (\\lambda w_1, \\lambda w_2, \\cdots, \\lambda w_n)\n",
    "$$\n",
    "将两边的各分量加在一起，并代入转移矩阵定义\n",
    "$$\n",
    "\\sum_{j = 1}^n p_{ij} = 1, \\quad i = 1, 2, \\cdots, n,\n",
    "$$\n",
    "得\n",
    "$$\n",
    "\\sum_{i = 1}^n w_i = \\lambda \\sum_{i = 1}^n w_i \\Rightarrow \\lambda = 1.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里注意，定理的要求是各元素均本质大于零，这是一个很强的条件。也就是任何两种状态，在一步之内，均有可能互变。实际上，这个条件可以放宽到：\n",
    "1. 系统的任何两个状态，都能在有限步互变，这一点又称为是**不可约**的；\n",
    "2. 系统的状态概率分布（注意不是状态）不存在周期变化；\n",
    "任何一个非周期不可约的马尔可夫链，都存在平稳分布。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "之前的推导已经明示了平稳分布的求法：\n",
    "+ 理论上，它可以通过计算转移矩阵的 $1$ 特征值对应的左特征向量得到；\n",
    "+ 计算上，它可以通过不断抽样模拟马尔可夫链（迭代）得到。 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC 抽样——Metropolis 算法\n",
    "\n",
    "我们继续抽样的话题。之前我们已经学习到：\n",
    "1. 如果存在累计分布函数的逆映射，则我们可以直接抽取指定分布；\n",
    "2. 如果只有概率密度函数，无法求得累计分布函数的逆映射，那么我们可以使用拒绝-接受法（AR）抽取指定分布；\n",
    "但是在一些实际工作中，我们只知道某个事件是随机的，并且知道它的一些分布规律，但是连概率密度函数 $f(x)$ 都无法给出。因为概率密度函数，必须满足\n",
    "$$\n",
    "\\int_{\\infty}^{\\infty} f(x) dx = 1,\n",
    "$$\n",
    "\n",
    "而很多时候，我们只知道某个随机事件，服从概率分布规律 $f(x)$，但是 $f(x)$ 无法写出显式表达式，或者有表达式，但是无法求积分，此时我们就没有办法对 $f(x)$ 做归一化，得到概率密度函数。或者说，在应用 AR 法时，本质困难出现在我们不知道 $f(x)$ 的上界，它甚至没有上界，就是一个一直增长的随机事件。\n",
    "\n",
    "在这种情况下，我们仍然可以利用马尔可夫链来获得指定随机事件的抽取（随机模拟）。基本想法是，我们构建一个马尔可夫链，使得它的平稳分布趋于指定的随机事件。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先，我们要假设，在指定的状态空间 $\\Omega = \\{x_1, x_2, \\cdots, x_n\\}$ 上，能够产生出我们的目标分布的概率密度函数，这里主要指状态空间要一致。然后关键在于，如果当前的状态是 $X_t$，那么下一个时刻的状态 $X_{t + 1}$ 如果确定，或者说，在已知当前状态为 $X_t$ 的前提下，$X_{t + 1}$ 满足什么概率分布。我们用概率论的记号，将这个分布记为 $g(\\cdot | X_t)$，称为**建议分布**。对应状态空间 $\\Omega$ 中的一个具体的状态 $x_i$, $i = 1, 2, \\cdots, n$，函数 $g$ 不必在整个 $\\Omega$ 上都本质大于零，实际上，$\\forall x \\in \\Omega$，我们称满足 $g(y | x) > 0$ 的状态构成的子集，为状态 $x$ 的**邻域**，记作 $\\mathcal{N}_x$。在很多算法中，我们要求 $g$ 满足对称性，也即\n",
    "$$\n",
    "g(y | x) = g(x | y), \\quad \\forall x, y \\in \\Omega.\n",
    "$$\n",
    "这一点在模型中往往会自然地做到（符合物理实际），并且它实际上可以弱化为\n",
    "$$\n",
    "g(y | x) > 0 \\Leftrightarrow g(x | y) > 0,  \\quad \\forall x, y \\in \\Omega,\n",
    "$$\n",
    "也即任两种状态的转移是可逆的。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去我们就可以逐一通过 $g$ 来模拟马尔可夫链的转移，但每一次从 $x$ 抽取到 $y$，我们都要做一个判定，到底需不需要转。也就是说，即便有 $X_t = x$，实际通过 $g(\\cdot | x)$ 抽取得到了 $y$，$X_{t + 1}$ 是否就取 $y$，还需要做一个独立抽取来判定，这里首先定义接受概率：\n",
    "$$\n",
    "h(x, y) = \\min \\left\\{1, \\frac{f(y)g(x | y)}{f(x)g(y | x)}\\right\\},\n",
    "$$\n",
    "特别地，当 $g$ 满足对称性条件时，上式简化为\n",
    "$$\n",
    "h(x, y) = \\min \\left\\{1, \\frac{f(y)}{f(x)}\\right\\}.\n",
    "$$\n",
    "我们首先计算这个概率，然后抽取一个 $(0, 1)$ 上的均匀分布的随机数 $r$，只有当 $r < h(x, y)$ 时，我们才令 $X_{t + 1} = y$，否则，仍然令 $X_{t + 1} = x$。也就是说，在这一个时间步，依概率 $h(x, y)$ 决定马尔可夫链是否从 $x$ 转移到 $y$。（这一个判定的意义在于，如果 $y$ 是一个概率极低的状态，那么它很难被转移到，但并不是完全没有希望。）\n",
    "这里，$f$ 只需要可计算就行，实际不需要满足概率密度函数归一化的要求（当然保正还是必须的）。而 $g$ 必须是一个严格的概率分布，但它是由我们设计的。这一个算法，被称为 Metropolis-Hastings 算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例** 我们先直接模拟两颗骰子的情形。显然两颗骰子的状态空间是 $\\Omega = \\{2, 3, \\cdots, 12\\}$。这个模型当然是可以得到概率质量分布的，但是我们这里直接通过 Metropolis-Hasting 算法直接产生符合这个模型的随机数。我们首先确定一下每一种状态可能的组合数，它实际上就对应这里的 $f(x)$，显然，\n",
    "$$\n",
    "f(2) = 1, f(3) = 2, f(4) = 3, f(5) = 4, f(6) = 5, f(7) = 6, f(8) = 5, f(9) = 4, f(10) = 3, f(11) = 2, f(12) = 1.\n",
    "$$\n",
    "注意这里 $f(x)$ 并未归一化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去考虑如何取邻域，也就是对 $x \\in \\Omega$，下一步 $y$ 可能的状态分布。这个取法可以很灵活，比如，取**极小邻域**，也就是 $\\forall x \\in \\Omega$，$\\mathcal{N}_x = \\{x - 1, x + 1\\}$，进一步加上对称性要求，我们就可以取\n",
    "$$\n",
    "g (y | x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\\frac{1}{2}, & y = \\max\\{x - 1, 2\\}, \\\\\n",
    "\\frac{1}{2}, & y = \\min\\{x + 1, 12\\},\\\\\n",
    "0, & \\mbox{others}.\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里注意考虑了 $\\Omega$ 两端的边界情况：\n",
    "$$\n",
    "g(3 | 2) = g(2 | 2) = g(11 / 12) = g(12 | 12) = \\frac{1}{2}. \n",
    "$$\n",
    "保证了对称性。这个自己设计的概率密度函数（严格说应该是概率质量分布）也可以等价地用下面的矩阵表示：\n",
    "$$\n",
    "G = \\left(\n",
    "\\begin{array}{ccccccc}\n",
    "\\frac{1}{2} & \\frac{1}{2} & 0 & \\cdots & 0 & 0 & 0 \\\\\n",
    "\\frac{1}{2} & 0 & \\frac{1}{2} & \\cdots & 0 & 0 & 0 \\\\\n",
    "0 & \\frac{1}{2} & 0 & \\cdots & 0 & 0 & 0 \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots & \\vdots \\\\\n",
    "0 & 0 & 0 & \\cdots & 0 & \\frac{1}{2} & 0 \\\\\n",
    "0 & 0 & 0 & \\cdots & \\frac{1}{2} & 0 & \\frac{1}{2} \\\\\n",
    "0 & 0 & 0 & \\cdots & 0 & \\frac{1}{2} & \\frac{1}{2}\n",
    "\\end{array}\n",
    "\\right).\n",
    "$$\n",
    "这里分布的对称性和矩阵的对称性是一致的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20000\n",
    "f = np.array([0, 0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1])  \n",
    "d = np.zeros(N)\n",
    "x = 5\n",
    "for i in range(N):\n",
    "    U = np.random.rand()\n",
    "    if x is 2:\n",
    "        if U < 0.5:\n",
    "            y = 3\n",
    "        else:\n",
    "            y = 2\n",
    "    elif x is 12:\n",
    "        if U < 0.5:\n",
    "            y = 11\n",
    "        else:\n",
    "            y = 12\n",
    "    else:\n",
    "        if U < 0.5:\n",
    "            y = x - 1\n",
    "        else:\n",
    "            y = x + 1\n",
    "    h = np.min([1, f[y] / f[x]])\n",
    "    U = np.random.rand()\n",
    "    if U < h:\n",
    "        x = y;\n",
    "    d[i] = x        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 618., 1184., 1706., 2214., 2746., 3241., 2729., 2227., 1633.,\n",
       "        1145.,  557.]),\n",
       " array([ 2.        ,  2.90909091,  3.81818182,  4.72727273,  5.63636364,\n",
       "         6.54545455,  7.45454545,  8.36363636,  9.27272727, 10.18181818,\n",
       "        11.09090909, 12.        ]),\n",
       " <a list of 11 Patch objects>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAEUdJREFUeJzt3X+s3XV9x/Hny4K/3YBxIdjWXeY6J5pZTINsJIsThQLGYjKWkk07R1L/gA0Xk624P3A6lpqpTDNlqdJRNwYSxNBIJ3bIYkzGj4IMKZVxBx1c29HrQMSR4Yrv/XG+jQe47T23vfccej/PR3Jyvt/3+XzP9/0J5b7u98c5N1WFJKk9Lxl1A5Kk0TAAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY06YtQNHMixxx5b4+Pjo25Dkg4rd9111w+qamymcS/qABgfH2fbtm2jbkOSDitJ/nOQcZ4CkqRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRr2oPwksDcP4upuGur+d688Z6v6k/fEIQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjZgyAJC9PckeSf0uyPcmfd/UTk9ye5MEkX07y0q7+sm59ont9vO+9LunqDyQ5c74mJUma2SBHAM8A76iqtwDLgZVJTgU+AVxeVcuAJ4ALuvEXAE9U1S8Dl3fjSHISsBp4E7AS+HySRXM5GUnS4GYMgOr5cbd6ZPco4B3A9V19E3But7yqW6d7/fQk6erXVtUzVfUwMAGcMiezkCTN2kDXAJIsSnIPsAfYCvwH8MOq2tsNmQQWd8uLgUcButefBH6hvz7NNpKkIRsoAKrq2apaDiyh91v7G6cb1j1nP6/tr/4cSdYm2ZZk29TU1CDtSZIOwqzuAqqqHwL/ApwKHJVk398TWALs6pYngaUA3es/DzzeX59mm/59bKiqFVW1YmxsbDbtSZJmYZC7gMaSHNUtvwJ4J7ADuBX47W7YGuDGbnlzt073+jerqrr66u4uoROBZcAdczURSdLsDPIXwU4ANnV37LwEuK6qvpbkfuDaJH8BfAe4sht/JfD3SSbo/ea/GqCqtie5Drgf2AtcWFXPzu10JEmDmjEAqupe4ORp6g8xzV08VfW/wHn7ea/LgMtm36Ykaa75SWBJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRg3wZnDRU4+tuGnULUhMMAGnIhhlwO9efM7R96fDjKSBJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjZoxAJIsTXJrkh1Jtie5uKt/NMn3k9zTPc7u2+aSJBNJHkhyZl99ZVebSLJufqYkSRrEIN8FtBf4cFXdneQ1wF1JtnavXV5Vn+wfnOQkYDXwJuC1wD8n+ZXu5c8B7wImgTuTbK6q++diIpKk2ZkxAKpqN7C7W34qyQ5g8QE2WQVcW1XPAA8nmQBO6V6bqKqHAJJc2401ACRpBGZ1DSDJOHAycHtXuijJvUk2Jjm6qy0GHu3bbLKr7a8uSRqBgQMgyauBrwAfqqofAVcArweW0ztC+NS+odNsXgeoP38/a5NsS7Jtampq0PYkSbM0UAAkOZLeD/+rq+oGgKp6rKqeraqfAl/gZ6d5JoGlfZsvAXYdoP4cVbWhqlZU1YqxsbHZzkeSNKBB7gIKcCWwo6o+3Vc/oW/Ye4H7uuXNwOokL0tyIrAMuAO4E1iW5MQkL6V3oXjz3ExDkjRbg9wFdBrwPuC7Se7pah8Bzk+ynN5pnJ3ABwGqanuS6+hd3N0LXFhVzwIkuQi4GVgEbKyq7XM4F0nSLAxyF9C3mf78/ZYDbHMZcNk09S0H2k6SNDx+EliSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRg3wdtMT4uptG3YIOwrD/u+1cf85Q96dD4xGAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1asYASLI0ya1JdiTZnuTirn5Mkq1JHuyej+7qSfLZJBNJ7k3y1r73WtONfzDJmvmbliRpJoMcAewFPlxVbwROBS5MchKwDrilqpYBt3TrAGcBy7rHWuAK6AUGcCnwNuAU4NJ9oSFJGr4ZA6CqdlfV3d3yU8AOYDGwCtjUDdsEnNstrwK+VD23AUclOQE4E9haVY9X1RPAVmDlnM5GkjSwWV0DSDIOnAzcDhxfVbuhFxLAcd2wxcCjfZtNdrX91SVJIzBwACR5NfAV4ENV9aMDDZ2mVgeoP38/a5NsS7Jtampq0PYkSbM0UAAkOZLeD/+rq+qGrvxYd2qH7nlPV58ElvZtvgTYdYD6c1TVhqpaUVUrxsbGZjMXSdIsDHIXUIArgR1V9em+lzYD++7kWQPc2Fd/f3c30KnAk90popuBM5Ic3V38PaOrSZJGYJA/CHMa8D7gu0nu6WofAdYD1yW5AHgEOK97bQtwNjABPA18AKCqHk/yceDObtzHqurxOZmFJGnWZgyAqvo205+/Bzh9mvEFXLif99oIbJxNg5Kk+eEngSWpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowb5ewB6ERpfd9OoW5B0mPMIQJIaZQBIUqMMAElqlNcAJM2ZYV+b2rn+nKHub6HxCECSGmUASFKjDABJapQBIEmNmjEAkmxMsifJfX21jyb5fpJ7usfZfa9dkmQiyQNJzuyrr+xqE0nWzf1UJEmzMcgRwFXAymnql1fV8u6xBSDJScBq4E3dNp9PsijJIuBzwFnAScD53VhJ0ojMeBtoVX0ryfiA77cKuLaqngEeTjIBnNK9NlFVDwEkubYbe/+sO5YkzYlDuQZwUZJ7u1NER3e1xcCjfWMmu9r+6pKkETnYALgCeD2wHNgNfKqrZ5qxdYD6CyRZm2Rbkm1TU1MH2Z4kaSYHFQBV9VhVPVtVPwW+wM9O80wCS/uGLgF2HaA+3XtvqKoVVbVibGzsYNqTJA3goAIgyQl9q+8F9t0htBlYneRlSU4ElgF3AHcCy5KcmOSl9C4Ubz74tiVJh2rGi8BJrgHeDhybZBK4FHh7kuX0TuPsBD4IUFXbk1xH7+LuXuDCqnq2e5+LgJuBRcDGqto+57ORJA1skLuAzp+mfOUBxl8GXDZNfQuwZVbdSZLmjZ8ElqRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEz/j0ADWZ83U2jbkGSZsUjAElqlAEgSY3yFJCkw9YwT73uXH/O0PY1LB4BSFKjDABJapQBIEmNMgAkqVEzBkCSjUn2JLmvr3ZMkq1JHuyej+7qSfLZJBNJ7k3y1r5t1nTjH0yyZn6mI0ka1CBHAFcBK59XWwfcUlXLgFu6dYCzgGXdYy1wBfQCA7gUeBtwCnDpvtCQJI3GjAFQVd8CHn9eeRWwqVveBJzbV/9S9dwGHJXkBOBMYGtVPV5VTwBbeWGoSJKG6GCvARxfVbsBuufjuvpi4NG+cZNdbX91SdKIzPVF4ExTqwPUX/gGydok25Jsm5qamtPmJEk/c7AB8Fh3aofueU9XnwSW9o1bAuw6QP0FqmpDVa2oqhVjY2MH2Z4kaSYHGwCbgX138qwBbuyrv7+7G+hU4MnuFNHNwBlJju4u/p7R1SRJIzLjdwEluQZ4O3Bskkl6d/OsB65LcgHwCHBeN3wLcDYwATwNfACgqh5P8nHgzm7cx6rq+ReWJUlDNGMAVNX5+3np9GnGFnDhft5nI7BxVt1JkuaNnwSWpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGjXj3wQ+nI2vu2nULUjSi5ZHAJLUKANAkhplAEhSow7pGkCSncBTwLPA3qpakeQY4MvAOLAT+J2qeiJJgM8AZwNPA79fVXcfyv4laViGfU1x5/pz5n0fc3EE8FtVtbyqVnTr64BbqmoZcEu3DnAWsKx7rAWumIN9S5IO0nycAloFbOqWNwHn9tW/VD23AUclOWEe9i9JGsChBkAB30hyV5K1Xe34qtoN0D0f19UXA4/2bTvZ1SRJI3ConwM4rap2JTkO2JrkewcYm2lq9YJBvSBZC/C6173uENuTJO3PIR0BVNWu7nkP8FXgFOCxfad2uuc93fBJYGnf5kuAXdO854aqWlFVK8bGxg6lPUnSARx0ACR5VZLX7FsGzgDuAzYDa7pha4Abu+XNwPvTcyrw5L5TRZKk4TuUU0DHA1/t3d3JEcA/VtXXk9wJXJfkAuAR4Lxu/BZ6t4BO0LsN9AOHsG9J0iE66ACoqoeAt0xT/2/g9GnqBVx4sPuTJM0tPwksSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqOGHgBJViZ5IMlEknXD3r8kqWeoAZBkEfA54CzgJOD8JCcNswdJUs+wjwBOASaq6qGq+glwLbBqyD1Ikhh+ACwGHu1bn+xqkqQhO2LI+8s0tXrOgGQtsLZb/XGSBw5hf8cCPziE7Q9Hrc25tfmCc25CPnFIc/7FQQYNOwAmgaV960uAXf0DqmoDsGEudpZkW1WtmIv3Oly0NufW5gvOuRXDmPOwTwHdCSxLcmKSlwKrgc1D7kGSxJCPAKpqb5KLgJuBRcDGqto+zB4kST3DPgVEVW0Btgxpd3NyKukw09qcW5svOOdWzPucU1Uzj5IkLTh+FYQkNWrBBUCSpUluTbIjyfYkF4+6p2FJsijJd5J8bdS9DEOSo5Jcn+R73X/vXx91T/MtyR93/67vS3JNkpePuqe5lmRjkj1J7uurHZNka5IHu+ejR9njXNvPnP+q+7d9b5KvJjlqrve74AIA2At8uKreCJwKXNjQ101cDOwYdRND9Bng61X1q8BbWOBzT7IY+CNgRVW9md6NFKtH29W8uApY+bzaOuCWqloG3NKtLyRX8cI5bwXeXFW/Bvw7cMlc73TBBUBV7a6qu7vlp+j9UFjwnzZOsgQ4B/jiqHsZhiQ/B/wmcCVAVf2kqn442q6G4gjgFUmOAF7J8z5HsxBU1beAx59XXgVs6pY3AecOtal5Nt2cq+obVbW3W72N3uem5tSCC4B+ScaBk4HbR9vJUPw18CfAT0fdyJD8EjAF/F132uuLSV416qbmU1V9H/gk8AiwG3iyqr4x2q6G5viq2g29X/KA40bcz7D9AfBPc/2mCzYAkrwa+Arwoar60aj7mU9J3g3sqaq7Rt3LEB0BvBW4oqpOBv6HhXda4Dm6896rgBOB1wKvSvJ7o+1K8y3Jn9E7tX31XL/3ggyAJEfS++F/dVXdMOp+huA04D1JdtL7htV3JPmH0bY07yaByarad3R3Pb1AWMjeCTxcVVNV9X/ADcBvjLinYXksyQkA3fOeEfczFEnWAO8Gfrfm4Z79BRcASULvvPCOqvr0qPsZhqq6pKqWVNU4vYuC36yqBf2bYVX9F/Bokjd0pdOB+0fY0jA8Apya5JXdv/PTWeAXvvtsBtZ0y2uAG0fYy1AkWQn8KfCeqnp6Pvax4AKA3m/D76P3W/A93ePsUTelefGHwNVJ7gWWA3854n7mVXe0cz1wN/Bdev//LrhPyCa5BvhX4A1JJpNcAKwH3pXkQeBd3fqCsZ85/w3wGmBr93Psb+d8v34SWJLatBCPACRJAzAAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElq1P8DLmuV2ERNKCIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(d, bins = 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "符合双骰子线性增长的规律。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "和极小邻域法类似，还有一种几乎万能的邻域取法是**极大邻域**法。这里我们直接取\n",
    "$$\n",
    "G = \\left(\n",
    "\\begin{array}{cccc}\n",
    "\\frac{1}{11} & \\frac{1}{11} & \\cdots & \\frac{1}{11} \\\\\n",
    "\\frac{1}{11} & \\frac{1}{11} & \\cdots & \\frac{1}{11} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\frac{1}{11} & \\frac{1}{11} & \\cdots & \\frac{1}{11} \\\\\n",
    "\\end{array}\n",
    "\\right),\n",
    "$$\n",
    "这里实际上取 $\\mathcal{N}_x = \\Omega$, 而 $g(\\cdot | x)$ 取了一个常量，这样的函数当然是对成的。这个模拟实际上写起来更简单："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 5\n",
    "for i in range(N):\n",
    "    y = np.random.randint(2, 13)\n",
    "    h = np.min([1, f[y] / f[x]])\n",
    "    U = np.random.rand()\n",
    "    if U < h:\n",
    "        x = y;\n",
    "    d[i] = x        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 545., 1158., 1665., 2285., 2876., 3160., 2856., 2148., 1585.,\n",
       "        1171.,  551.]),\n",
       " array([ 2.        ,  2.90909091,  3.81818182,  4.72727273,  5.63636364,\n",
       "         6.54545455,  7.45454545,  8.36363636,  9.27272727, 10.18181818,\n",
       "        11.09090909, 12.        ]),\n",
       " <a list of 11 Patch objects>)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAETtJREFUeJzt3W2MXddd7/HvD6ct0BaSkkkUbJcJYKApokllpeFWuioNTZwE4VSikiNorRLJvEigRZXAgRfhKVdGF9pLRQkKxNSF0BD1QbEa09SEogqJtnFKSOOYkCE1ydQmnpI2LVS3XIf/fXGWxakznjljnzknmfX9SEdn7/9ee++15PH8Zj+dk6pCktSfb5l2ByRJ02EASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjp11rQ7sJRzzz23Zmdnp90NSXpBeeCBB75UVTPLtXteB8Ds7CwHDhyYdjck6QUlyb+M0s5TQJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1Knn9ZPA0iTM7rxnovs7vOuaie5POhWPACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnlg2AJN+a5LNJ/iHJwSS/3uoXJvlMkseS/EWSF7f6S9r8XFs+O7Stm1r90SRXrtagJEnLG+UI4BvAG6vqNcDFwJYklwG/DbynqjYBXwaub+2vB75cVd8PvKe1I8lFwDbg1cAW4A+SrBvnYCRJo1s2AGrg39vsi9qrgDcCH2r1PcC1bXprm6ctvzxJWv3OqvpGVX0BmAMuHcsoJEkrNtI1gCTrkjwIHAP2A/8MfKWqjrcm88D6Nr0eeBKgLX8G+K7h+iLrDO9rR5IDSQ4sLCysfESSpJGMFABV9WxVXQxsYPBX+6sWa9bec4plp6qfvK/bqmpzVW2emVn2S+0lSadpRXcBVdVXgL8BLgPOTnLis4Q2AEfa9DywEaAt/07g6eH6IutIkiZslLuAZpKc3aa/Dfhx4BDwSeCnWrPtwN1tem+bpy3/66qqVt/W7hK6ENgEfHZcA5EkrcwonwZ6AbCn3bHzLcBdVfWxJI8Adyb5LeDvgdtb+9uBP00yx+Av/20AVXUwyV3AI8Bx4Iaqena8w5EkjWrZAKiqh4BLFqk/ziJ38VTV/wXecopt3QLcsvJuSpLGzSeBJalTBoAkdcpvBNPzzqS/oUvqlUcAktQpjwCkCZvkEY7fP6yleAQgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTywZAko1JPpnkUJKDSd7R6r+W5ItJHmyvq4fWuSnJXJJHk1w5VN/SanNJdq7OkCRJoxjlO4GPA++qqs8leTnwQJL9bdl7qup3hhsnuQjYBrwa+G7gr5L8QFv8PuBNwDxwf5K9VfXIOAYiSVqZZQOgqo4CR9v015IcAtYvscpW4M6q+gbwhSRzwKVt2VxVPQ6Q5M7W1gCQpClY0TWAJLPAJcBnWunGJA8l2Z3knFZbDzw5tNp8q52qLkmagpEDIMnLgA8D76yqrwK3At8HXMzgCOF3TzRdZPVaon7yfnYkOZDkwMLCwqjdkySt0EgBkORFDH7531FVHwGoqqeq6tmq+i/gj/jv0zzzwMah1TcAR5aof5Oquq2qNlfV5pmZmZWOR5I0olHuAgpwO3Coqt49VL9gqNmbgYfb9F5gW5KXJLkQ2AR8Frgf2JTkwiQvZnCheO94hiFJWqlR7gJ6PfBW4PNJHmy1XwGuS3Ixg9M4h4GfA6iqg0nuYnBx9zhwQ1U9C5DkRuBeYB2wu6oOjnEskqQVGOUuoL9l8fP3+5ZY5xbglkXq+5ZaT5I0OT4JLEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnRrloyAkZnfeM+0uSBozjwAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnfJBMGkNm/QDfId3XTPR/enMeAQgSZ0yACSpUwaAJHXKAJCkThkAktSpZQMgycYkn0xyKMnBJO9o9Vck2Z/ksfZ+TqsnyXuTzCV5KMlrh7a1vbV/LMn21RuWJGk5oxwBHAfeVVWvAi4DbkhyEbATuK+qNgH3tXmAq4BN7bUDuBUGgQHcDLwOuBS4+URoSJImb9kAqKqjVfW5Nv014BCwHtgK7GnN9gDXtumtwAdq4NPA2UkuAK4E9lfV01X1ZWA/sGWso5EkjWxF1wCSzAKXAJ8Bzq+qozAICeC81mw98OTQavOtdqq6JGkKRg6AJC8DPgy8s6q+ulTTRWq1RP3k/exIciDJgYWFhVG7J0laoZECIMmLGPzyv6OqPtLKT7VTO7T3Y60+D2wcWn0DcGSJ+jepqtuqanNVbZ6ZmVnJWCRJKzDKXUABbgcOVdW7hxbtBU7cybMduHuo/rZ2N9BlwDPtFNG9wBVJzmkXf69oNUnSFIzyYXCvB94KfD7Jg632K8Au4K4k1wNPAG9py/YBVwNzwNeBtwNU1dNJfhO4v7X7jap6eiyjkCSt2LIBUFV/y+Ln7wEuX6R9ATecYlu7gd0r6aAkaXX4JLAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1apSvhNTz0OzOe6bdBUkvcB4BSFKnDABJ6pQBIEmdMgAkqVNeBJY0NpO+OeHwrmsmur+1xiMASerUsgGQZHeSY0keHqr9WpIvJnmwva4eWnZTkrkkjya5cqi+pdXmkuwc/1AkSSsxyhHA+4Eti9TfU1UXt9c+gCQXAduAV7d1/iDJuiTrgPcBVwEXAde1tpKkKVn2GkBVfSrJ7Ijb2wrcWVXfAL6QZA64tC2bq6rHAZLc2do+suIeS5LG4kyuAdyY5KF2iuicVlsPPDnUZr7VTlV/jiQ7khxIcmBhYeEMuidJWsrpBsCtwPcBFwNHgd9t9SzStpaoP7dYdVtVba6qzTMzM6fZPUnSck7rNtCqeurEdJI/Aj7WZueBjUNNNwBH2vSp6pKkKTitI4AkFwzNvhk4cYfQXmBbkpckuRDYBHwWuB/YlOTCJC9mcKF47+l3W5J0ppY9AkjyQeANwLlJ5oGbgTckuZjBaZzDwM8BVNXBJHcxuLh7HLihqp5t27kRuBdYB+yuqoNjH40kaWSj3AV03SLl25dofwtwyyL1fcC+FfVOkrRqfBJYkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktSp0/pOYEl6Ppjdec/E9nV41zUT29ekGABjMskfREkaB08BSVKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4tGwBJdic5luThodorkuxP8lh7P6fVk+S9SeaSPJTktUPrbG/tH0uyfXWGI0ka1ShHAO8HtpxU2wncV1WbgPvaPMBVwKb22gHcCoPAAG4GXgdcCtx8IjQkSdOxbABU1aeAp08qbwX2tOk9wLVD9Q/UwKeBs5NcAFwJ7K+qp6vqy8B+nhsqkqQJOt1rAOdX1VGA9n5eq68HnhxqN99qp6o/R5IdSQ4kObCwsHCa3ZMkLWfcF4GzSK2WqD+3WHVbVW2uqs0zMzNj7Zwk6b+dbgA81U7t0N6Ptfo8sHGo3QbgyBJ1SdKUnG4A7AVO3MmzHbh7qP62djfQZcAz7RTRvcAVSc5pF3+vaDVJ0pQs+2mgST4IvAE4N8k8g7t5dgF3JbkeeAJ4S2u+D7gamAO+DrwdoKqeTvKbwP2t3W9U1ckXliVJE7RsAFTVdadYdPkibQu44RTb2Q3sXlHvJEmrxieBJalTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASerUGQVAksNJPp/kwSQHWu0VSfYneay9n9PqSfLeJHNJHkry2nEMQJJ0es4awzZ+rKq+NDS/E7ivqnYl2dnmfxm4CtjUXq8Dbm3vq2Z25z2ruXlJHZn075PDu65Z9X2sximgrcCeNr0HuHao/oEa+DRwdpILVmH/kqQRnGkAFPCJJA8k2dFq51fVUYD2fl6rrweeHFp3vtUkSVNwpqeAXl9VR5KcB+xP8o9LtM0itXpOo0GQ7AB45StfeYbdkySdyhkdAVTVkfZ+DPgocCnw1IlTO+39WGs+D2wcWn0DcGSRbd5WVZuravPMzMyZdE+StITTDoAkL03y8hPTwBXAw8BeYHtrth24u03vBd7W7ga6DHjmxKkiSdLknckpoPOBjyY5sZ0/r6qPJ7kfuCvJ9cATwFta+33A1cAc8HXg7Wewb0nSGTrtAKiqx4HXLFL/N+DyReoF3HC6+5MkjZdPAktSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjo18QBIsiXJo0nmkuyc9P4lSQMTDYAk64D3AVcBFwHXJblokn2QJA1M+gjgUmCuqh6vqv8E7gS2TrgPkiQmHwDrgSeH5udbTZI0YWdNeH9ZpFbf1CDZAexos/+e5NEz2N+5wJfOYP0Xot7G3Nt4wTF3Ib99RmP+nlEaTToA5oGNQ/MbgCPDDarqNuC2cewsyYGq2jyObb1Q9Dbm3sYLjrkXkxjzpE8B3Q9sSnJhkhcD24C9E+6DJIkJHwFU1fEkNwL3AuuA3VV1cJJ9kCQNTPoUEFW1D9g3od2N5VTSC0xvY+5tvOCYe7HqY05VLd9KkrTm+FEQktSpNRcASTYm+WSSQ0kOJnnHtPs0KUnWJfn7JB+bdl8mIcnZST6U5B/bv/ePTrtPqy3JL7af64eTfDDJt067T+OWZHeSY0keHqq9Isn+JI+193Om2cdxO8WY/3f72X4oyUeTnD3u/a65AACOA++qqlcBlwE3dPRxE+8ADk27ExP0e8DHq+qHgNewxseeZD3wC8DmqvphBjdSbJtur1bF+4EtJ9V2AvdV1Sbgvja/lryf5455P/DDVfUjwD8BN417p2suAKrqaFV9rk1/jcEvhTX/tHGSDcA1wB9Puy+TkOQ7gP8J3A5QVf9ZVV+Zbq8m4izg25KcBXw7Jz1HsxZU1aeAp08qbwX2tOk9wLUT7dQqW2zMVfWJqjreZj/N4LmpsVpzATAsySxwCfCZ6fZkIv4P8EvAf027IxPyvcAC8CfttNcfJ3nptDu1mqrqi8DvAE8AR4FnquoT0+3VxJxfVUdh8EcecN6U+zNpPwv85bg3umYDIMnLgA8D76yqr067P6spyU8Ax6rqgWn3ZYLOAl4L3FpVlwD/wdo7LfBN2nnvrcCFwHcDL03yM9PtlVZbkl9lcGr7jnFve00GQJIXMfjlf0dVfWTa/ZmA1wM/meQwg09YfWOSP5tul1bdPDBfVSeO7j7EIBDWsh8HvlBVC1X1/4CPAP9jyn2alKeSXADQ3o9NuT8TkWQ78BPAT9cq3LO/5gIgSRicFz5UVe+edn8moapuqqoNVTXL4KLgX1fVmv7LsKr+FXgyyQ+20uXAI1Ps0iQ8AVyW5Nvbz/nlrPEL30P2Atvb9Hbg7in2ZSKSbAF+GfjJqvr6auxjzQUAg7+G38rgr+AH2+vqaXdKq+LngTuSPARcDPyvKfdnVbWjnQ8BnwM+z+D/75p7QjbJB4G/A34wyXyS64FdwJuSPAa8qc2vGacY8+8DLwf2t99jfzj2/foksCT1aS0eAUiSRmAASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUqf8PgqyBDa1CEHAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(d, bins = 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个算法，也可以使用在连续分布抽样上。比如，已知密度函数为\n",
    "$$\n",
    "f(x) = 0.5 x^2 e^{-x}\n",
    "$$\n",
    "的伽马分布，这里我们不用考察这个密度函数是否归一，它的状态空间是 $\\Omega = [0, +\\infty)$，所以我们在产生建议分布 $g$ 时也要注意，$g(y | x)$ 必须确保 $y > 0$，同时我们也希望继续保持 $g$ 的对称性。为此，我们取 $\\mathcal{N}_x = (x - 1, x + 1)$, 且 $g(\\cdot | x)$ 为其上的均匀分布。如果这里 $y < 0$，则令 $y = x$，这个做法类似之前离散的情况。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 0.5 * x * x * np.exp(-x) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 40000\n",
    "d = np.zeros(N)\n",
    "x = 2\n",
    "for i in range(N):\n",
    "    y = np.random.rand() * 2 - 1 + x\n",
    "    if y < 0:\n",
    "        y = x\n",
    "    h = np.min([1, f(y) / f(x)])\n",
    "    U = np.random.rand()\n",
    "    if U < h:\n",
    "        x = y\n",
    "    d[i] = x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD8CAYAAACW/ATfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAEsFJREFUeJzt3X+s3Xd93/HnawkphUKdkJvUs52arlbaqlp+7ArSRao6XDYSUJw/miqsA5N58v4IDNZOjVml/ZCmydWmpiCmTBahOF0GZCkoFssolgFVk5aM6xDCD8PiMhrf2o3dEkLXaGMZ7/1xPi4H+zr3e3/53HM/z4d09P1+P+dzznlf6/p1PudzPt/vTVUhSdrY/sqkC5AkrT3DXpI6YNhLUgcMe0nqgGEvSR0w7CWpA4a9JHXAsJekDhj2ktSBSyddAMCVV15Z27dvn3QZkjRVjh49+qdVNTOk77oI++3btzM3NzfpMiRpqiT5o6F9ncaRpA4Y9pLUAcNekjpg2EtSBwx7SeqAYS9JHTDsJakDhr0kdcCwl6QOrIszaDeC7fv+81/uf3P/mydYiSSdz5G9JHVg0bBPcm2SJ8du30nyniRXJDmc5Om2vbz1T5L3Jzme5KkkN679jyFJeimLhn1Vfb2qrq+q64G/AbwAfALYBxypqh3AkXYMcAuwo932AvetReGSpOGWOo2zE/jDqvojYBdwsLUfBG5v+7uAB2rkMWBTks2rUq0kaVmWGvZ3Ah9p+1dX1SmAtr2qtW8BTow9Zr61/YAke5PMJZk7c+bMEsuQJC3F4LBPchlwG/CfFuu6QFud11B1oKpmq2p2ZmbQtfclScu0lKWXtwBPVNWz7fjZJJur6lSbpjnd2ueBbWOP2wqcXHmp08NlmJLWm6WE/Vv5/hQOwCFgN7C/bR8Za39nko8CrweePzvdo+/zDUHSxTQo7JO8Angj8A/HmvcDDyXZAzwD3NHaHwVuBY4zWrlz16pVK0lalkFhX1UvAK85p+3PGK3OObdvAXevSnWSpFXhGbSS1AGvjXMRjc/TS9LF5Mhekjpg2EtSB5zGWQdchilprRn2S2QwS5pGTuNIUgcMe0nqgNM465hTRpJWiyN7SeqAYS9JHTDsJakDztmvsaVeIsFLKkhaC47sJakDhr0kdcBpnBVwykXStHBkL0kdMOwlqQOGvSR1wLCXpA4MCvskm5I8nORrSY4l+bkkVyQ5nOTptr289U2S9yc5nuSpJDeu7Y8gSVrM0JH9+4BPVdVPAdcBx4B9wJGq2gEcaccAtwA72m0vcN+qVixJWrJFl14meTXw88A7AKrqu8B3k+wCfqF1Owh8DrgH2AU8UFUFPNY+FWyuqlOrXn2nvBqmpKUass7+J4AzwO8kuQ44CrwbuPpsgFfVqSRXtf5bgBNjj59vbT8Q9kn2Mhr5c80116zkZ1hzrqeXNO2GTONcCtwI3FdVNwB/wfenbBaSBdrqvIaqA1U1W1WzMzMzg4qVJC3PkLCfB+ar6vF2/DCj8H82yWaAtj091n/b2OO3AidXp1xJ0nIsGvZV9SfAiSTXtqadwFeBQ8Du1rYbeKTtHwLe3lbl3AQ873y9JE3W0GvjvAt4MMllwDeAuxi9UTyUZA/wDHBH6/socCtwHHih9ZUkTdCgsK+qJ4HZBe7auUDfAu5eYV2SpFXkGbSS1AHDXpI6YNhLUgcMe0nqgH+pakp4Fq+klXBkL0kdMOwlqQOGvSR1wLCXpA74Be2U89r2koZwZC9JHTDsJakDhr0kdcCwl6QOGPaS1AFX42wgrsyRdCGO7CWpA4a9JHXAaZwxG+nKkk7pSBo3aGSf5JtJvpTkySRzre2KJIeTPN22l7f2JHl/kuNJnkpy41r+AJKkxS1lGudvVdX1VXX2D4/vA45U1Q7gSDsGuAXY0W57gftWq1hJ0vKsZM5+F3Cw7R8Ebh9rf6BGHgM2Jdm8gteRJK3Q0LAv4NNJjibZ29qurqpTAG17VWvfApwYe+x8a5MkTcjQL2hvrqqTSa4CDif52kv0zQJtdV6n0ZvGXoBrrrlmYBmSpOUYNLKvqpNtexr4BPA64Nmz0zNte7p1nwe2jT18K3Bygec8UFWzVTU7MzOz/J9AkrSoRcM+ySuTvOrsPvC3gS8Dh4Ddrdtu4JG2fwh4e1uVcxPw/NnpHknSZAyZxrka+ESSs/3/Y1V9KsnngYeS7AGeAe5o/R8FbgWOAy8Ad6161ZKkJVk07KvqG8B1C7T/GbBzgfYC7l6V6iRJq8IzaDvg2bSSDHsB518qwjcFaWPxQmiS1AHDXpI6YNhLUgcMe0nqgGEvSR0w7CWpAy697Ixr7qU+ObKXpA4Y9pLUAcNekjpg2EtSBwx7SeqAYS9JHTDsJakDrrPv2LmXNZa0cTmyl6QOGPaS1AHDXpI6MDjsk1yS5AtJPtmOX5vk8SRPJ/lYksta+w+14+Pt/u1rU7okaailjOzfDRwbO/5N4N6q2gE8B+xp7XuA56rqJ4F7Wz9J0gQNCvskW4E3Ax9sxwHeADzcuhwEbm/7u9ox7f6drb8kaUKGjux/G/h14Hvt+DXAt6vqxXY8D2xp+1uAEwDt/udbf0nShCwa9kneApyuqqPjzQt0rQH3jT/v3iRzSebOnDkzqFhJ0vIMOanqZuC2JLcCLwdezWikvynJpW30vhU42frPA9uA+SSXAj8KfOvcJ62qA8ABgNnZ2fPeDDRZ/pETaWNZdGRfVe+tqq1VtR24E/hMVf0K8Fngl1q33cAjbf9QO6bd/5mqMswlaYJWss7+HuBXkxxnNCd/f2u/H3hNa/9VYN/KSpQkrVTWw6B7dna25ubmJvLaXh9maZzSkdaPJEeranZIXy+EpmVzXl+aHl4uQZI6YNhLUgcMe0nqgGEvSR3wC1otiauXpOnkyF6SOmDYS1IHDHtJ6oBz9lp1nmwlrT+O7CWpA4a9JHXAsJekDhj2ktQBw16SOtDdahzPAJXUI0f2ktQBw16SOmDYS1IHDHtJ6sCiYZ/k5Un+e5IvJvlKkn/Z2l+b5PEkTyf5WJLLWvsPtePj7f7ta/sjSJIWM2Rk/3+AN1TVdcD1wJuS3AT8JnBvVe0AngP2tP57gOeq6ieBe1s/SdIELbr0sqoK+F/t8GXtVsAbgL/b2g8C/wK4D9jV9gEeBj6QJO15tEG5pFVa3wbN2Se5JMmTwGngMPCHwLer6sXWZR7Y0va3ACcA2v3PA69ZzaIlSUszKOyr6v9V1fXAVuB1wE8v1K1t8xL3/aUke5PMJZk7c+bM0HolScuwpNU4VfVt4HPATcCmJGengbYCJ9v+PLANoN3/o8C3FniuA1U1W1WzMzMzy6tekjTIkNU4M0k2tf0fBn4ROAZ8Fvil1m038EjbP9SOafd/xvl6SZqsIdfG2QwcTHIJozeHh6rqk0m+Cnw0yb8CvgDc3/rfD/xukuOMRvR3rkHdkqQlGLIa5ynghgXav8Fo/v7c9v8N3LEq1UmSVkV3V73UxeXfo5XWB8NeF43BL02O18aRpA44stdEOMqXLi5H9pLUAcNekjpg2EtSBwx7SeqAYS9JHTDsJakDhr0kdcCwl6QOeFKVJs4TrKS158hekjpg2EtSBwx7SeqAYS9JHTDsJakDhr0kdcCwl6QOLLrOPsk24AHgx4DvAQeq6n1JrgA+BmwHvgn8clU9lyTA+4BbgReAd1TVE2tTvjYa19xLa2PISVUvAr9WVU8keRVwNMlh4B3Akaran2QfsA+4B7gF2NFurwfua9uJGQ8QSerRotM4VXXq7Mi8qv4cOAZsAXYBB1u3g8DtbX8X8ECNPAZsSrJ51SuXJA22pDn7JNuBG4DHgaur6hSM3hCAq1q3LcCJsYfNt7Zzn2tvkrkkc2fOnFl65ZKkwQaHfZIfAX4PeE9Vfeelui7QVuc1VB2oqtmqmp2ZmRlahiRpGQaFfZKXMQr6B6vq46352bPTM217urXPA9vGHr4VOLk65UqSlmPRsG+ra+4HjlXVb43ddQjY3fZ3A4+Mtb89IzcBz5+d7pEkTcaQ1Tg3A28DvpTkydb2T4H9wENJ9gDPAHe0+x5ltOzyOKOll3etasXqhsswpdWzaNhX1X9l4Xl4gJ0L9C/g7hXWJf0Ag19aGc+glaQO+JeqNHUc5UtL58hekjpg2EtSBzbsNI7Xw5Gk73NkL0kdMOwlqQMbdhpHfXPFjvSDHNlLUgcMe0nqgNM4mmpO10jDGPbaMFxuK12Y0ziS1AHDXpI6YNhLUgcMe0nqgGEvSR0w7CWpAy69VFdcl69eLTqyT/KhJKeTfHms7Yokh5M83baXt/YkeX+S40meSnLjWhYvSRpmyMj+w8AHgAfG2vYBR6pqf5J97fge4BZgR7u9HrivbaWJ8WQraUDYV9UfJNl+TvMu4Bfa/kHgc4zCfhfwQFUV8FiSTUk2V9Wp1SpYWmtO9WgjWu6c/dVnA7yqTiW5qrVvAU6M9ZtvbYa91h1DXT1Z7dU4WaCtFuyY7E0yl2TuzJkzq1yGJGnccsP+2SSbAdr2dGufB7aN9dsKnFzoCarqQFXNVtXszMzMMsuQJA2x3GmcQ8BuYH/bPjLW/s4kH2X0xezzztdrmp375a7TPZpWi4Z9ko8w+jL2yiTzwD9nFPIPJdkDPAPc0bo/CtwKHAdeAO5ag5qlVeeKHW10Q1bjvPUCd+1coG8Bd6+0KGm9utCXun7Zq/XOyyVIUgcMe0nqgNfGkZbJeX5NE0f2ktQBw16SOuA0jrSGXKWj9cKRvSR1YEON7P3CTJIWtqHCXloPLjTocEpHk2TYSxNg8Otic85ekjrgyF5aRxzxa60Y9tIU8E1AK2XYS+uUq8u0mgx7acIMdV0Mhr00ZS705uD19fVSXI0jSR2Y+pG9H4Gll+YoX+DIXpK6MPUje0kjQz7lOsrv15qM7JO8KcnXkxxPsm8tXkOSNNyqj+yTXAL8O+CNwDzw+SSHquqrq/1akpZvyCh/yMofTYe1mMZ5HXC8qr4BkOSjwC7AsJfWKRc6bHxrEfZbgBNjx/PA69fgdSStE0td+z+kz5BPD0v9dLKS55z27ztSVav7hMkdwN+pqn/Qjt8GvK6q3nVOv73A3nZ4LfD1ZbzclcCfrqDcSZnGuqexZpjOuqexZrDui+lszT9eVTNDHrAWI/t5YNvY8Vbg5LmdquoAcGAlL5RkrqpmV/IckzCNdU9jzTCddU9jzWDdF9Nyal6L1TifB3YkeW2Sy4A7gUNr8DqSpIFWfWRfVS8meSfw+8AlwIeq6iur/TqSpOHW5KSqqnoUeHQtnvscK5oGmqBprHsaa4bprHsaawbrvpiWXPOqf0ErSVp/vDaOJHVgasN+2i7JkGRbks8mOZbkK0nePemaliLJJUm+kOSTk65liCSbkjyc5Gvt3/znJl3TEEn+cfv9+HKSjyR5+aRrWkiSDyU5neTLY21XJDmc5Om2vXySNZ7rAjX/m/Y78lSSTyTZNMkaF7JQ3WP3/ZMkleTKxZ5nKsN+7JIMtwA/A7w1yc9MtqpFvQj8WlX9NHATcPcU1Dzu3cCxSRexBO8DPlVVPwVcxxTUnmQL8I+A2ar6WUYLHO6cbFUX9GHgTee07QOOVNUO4Eg7Xk8+zPk1HwZ+tqr+OvA/gPde7KIG+DDn102SbYwuS/PMkCeZyrBn7JIMVfVd4OwlGdatqjpVVU+0/T9nFD5bJlvVMEm2Am8GPjjpWoZI8mrg54H7Aarqu1X17clWNdilwA8nuRR4BQuco7IeVNUfAN86p3kXcLDtHwRuv6hFLWKhmqvq01X1Yjt8jNF5QevKBf6tAe4Ffh0Y9MXrtIb9QpdkmIrgBEiyHbgBeHyylQz224x+qb436UIG+gngDPA7berpg0leOemiFlNVfwz8W0YjtVPA81X16clWtSRXV9UpGA1ugKsmXM9S/X3gv0y6iCGS3Ab8cVV9cehjpjXss0DbVCwrSvIjwO8B76mq70y6nsUkeQtwuqqOTrqWJbgUuBG4r6puAP6C9TelcJ42x70LeC3wV4FXJvl7k62qD0l+g9FU64OTrmUxSV4B/Abwz5byuGkN+0GXZFhvkryMUdA/WFUfn3Q9A90M3Jbkm4ymy96Q5D9MtqRFzQPzVXX2k9PDjMJ/vftF4H9W1Zmq+r/Ax4G/OeGaluLZJJsB2vb0hOsZJMlu4C3Ar9R0rEX/a4wGBF9s/y+3Ak8k+bGXetC0hv3UXZIhSRjNIR+rqt+adD1DVdV7q2prVW1n9O/8mapa16PNqvoT4ESSa1vTTqbjEtvPADcleUX7fdnJFHyxPOYQsLvt7wYemWAtgyR5E3APcFtVvTDpeoaoqi9V1VVVtb39v5wHbmy/9xc0lWHfvlA5e0mGY8BDU3BJhpuBtzEaGT/ZbrdOuqgN7F3Ag0meAq4H/vWE61lU+yTyMPAE8CVG/z/X5dmdST4C/Dfg2iTzSfYA+4E3Jnma0SqR/ZOs8VwXqPkDwKuAw+3/5L+faJELuEDdS3+e6fjUIklaiakc2UuSlsawl6QOGPaS1AHDXpI6YNhLUgcMe0nqgGEvSR0w7CWpA/8fkCprQYvsoBgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "num = plt.hist(d[20000:], bins = 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对比一下目标曲线，下图中的蓝线是根据统计结果做的归一化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f95f9de8550>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNX9//HXZ2ayJ2QjYUkCAQXZBQmI4tIWF9QKaq3iUnGptrVWrV2+tlptbevPtt8u9qutWutSqyIuVaooouIKaMK+Q9hDgAQSQsg+M+f3x5lAiAEmZLmzfJ6PxzzuZObemQ+azHvuOeeeI8YYlFJKRR+X0wUopZRyhgaAUkpFKQ0ApZSKUhoASikVpTQAlFIqSmkAKKVUlNIAUEqpKKUBoJRSUUoDQCmlopTH6QJa69mzp8nPz3e6DKWUCiuLFi3aY4zJas8xIRcA+fn5FBUVOV2GUkqFFRHZ2t5jtAlIKaWilAaAUkpFKQ0ApZSKUhoASikVpTQAlFIqSmkAKKVUlNIAUEqpKBVy1wFEhYrNsO5tSM6GvPGQmgciTlellIoyGgDdxe+Hxc/A0hegpPDw51L6woTvwWm3gUtPypRS3UMDoDv4/fDWD2HRM5A9HM75JQy/FOr2wfYvYN1smPsL2DQPLn3cnhkopVQX0wDoasbA7B/ZD/8zfgiT7j/U3JMO9B0N42+GRU/DOz+Dv0+EK56F/qc7WbVSKgpoe0NXMgZm/xiKnoKJdx7+4d+SCBTcCDfPg/hUeOFKKFvT/fUqpaKKBkBXWvEKFD4Jp99um32O1dHbaxhc9zrEJMLzV0D17u6oUikVpTQAukpjDcy9D/qMhnN+Ffwon9RcuHoG1O6BF6dBY23X1qmUiloaAF3l079AdSlc8Lv2j+zpOwa+8U8oXQL/vaNr6lNKRT0NgK6wbxvM/yuMuBz6TTi+1xhyIXzlblgxE9a+1bn1KaUUGgBdY+59gMC5v+rY65z5I+g1Et68C+oqO6U0pZRqpgHQ2XYsglX/gYl32Pb8jnDHwNRHoKYc5tzTOfUppVSABkBnK3wKYpPh9Ns65/X6joYz7oSlz8OG9zrnNZVSCg2AzlVfBategxHfgLiUznvds34KPU+Ct+4Cb0Pnva5SKqoFFQAiMllE1olIsYjc3cbzd4nIahFZLiLvi0j/Fs/5RGRp4DarM4sPOStehqZaGDu9c183Jh4mPwj7ttrrCpRSqhMcMwBExA08ClwADAOuEpFhrXZbAhQYY0YBrwC/b/FcnTFmdOA2pZPqDk2LnoXeI6HvKZ3/2ieeAwO/Ch/9XjuElVKdIpgzgPFAsTFmkzGmEZgBTG25gzFmnjGm+YqlhUAHez/DUOkS2LUcTpnedVM7n/dr28z0yR+75vWVUlElmADIAba3+Lkk8NiR3AS83eLneBEpEpGFInJJWweIyC2BfYrKy8uDKCkELXoGPAkw6oque4/eI+Hkq+Dzx6Fya9e9j1IqKgQTAG19nTVt7ihyLVAA/KHFw/2MMQXA1cBfROSEL72YMU8YYwqMMQVZWVlBlBRiGg7YeX9GXGYnc+tKX7sXxAUf/KZr30cpFfGCCYASIK/Fz7lAaeudROQc4B5gijHm4FAVY0xpYLsJ+BAY04F6Q9Pat6DxAJxyXde/V2oOnPod2+Fcvr7r308pFbGCCYBCYJCIDBCRWGAacNhoHhEZAzyO/fAva/F4uojEBe73BCYCqzur+JCx/m1I7g2547vn/U77AXjitS9AKdUhxwwAY4wXuA2YA6wBZhpjVonIAyLSPKrnD0Ay8HKr4Z5DgSIRWQbMAx4yxkRWAPiaoPh9GHxe9y3nmJxl1w9Y8TJUbOqe91RKRZygVgQzxswGZrd67L4W9885wnHzgZEdKTDkbVsADfth8OTufd+Jt9trAj75k50uQiml2kmvBO6o9XPAHQcDzu7e903pbfsclr1oZx9VSql20gDoqPXvwIAzIS65+9/7jDsBsWsPKKVUO2kAdMSeYthb3P3NP81Sc2H01bDk33AgTK+fUEo5RgOgI9a/Y7eDznOuhtNuA1+DzhGklGo3DYCOWP8OZA+D9P7H3rerZA22ZyCF/4CmOufqUEqFHQ2A41W3z44AGny+05XYs4DavbBshtOVKKXCiAbA8do0D/xe59r/W8o/A/qcDAseBb/f6WqUUmFCA+B4bfnUrvyVU+B0JXb20dN+AHs3QPFcp6tRSoUJDYDjtW0h5I0Hd1DX0nW94ZdAjxyY/39OV6KUChMaAMejbh/sXgX9TnO6kkPcMXaSuC2fwK6VTlejlAoDGgDHY/sXgAmtAAAY8y07SVzhP5yuRCkVBjQAjse2+eCKgZyxTldyuMQMGPlNWD5Tl41USh2TBsDx2LYQ+o6G2ESnK/my8bfYhemXPO90JUqpEKcB0F5N9bBjEfSb4HQlbeszyjZNFf4D/D6nq1FKhTANgPYqXQK+xtBr/29p/C1QuQWK33O6EqVUCNMAaK9t8+02L0TPAACGXgwpfezi8UopdQQaAO21bSH0PAmSMp2u5MjcMXbFsI3vw96NTlejlApRGgDt4ffDts+hfwg3/zQ75TpweWDR005XopQKURoA7VG2GhqqQrv9v1lKbxhykR0N1FTvdDVKqRCkAdAe2xbYbTgEANhmoLoKWP2G05UopUKQBkB7lC6BpCxI6+d0JcEZcDZknghF/3S6EqVUCNIAaI+dy6DPaDv7ZjgQgbE3wPbPdX4gpdSXaAAEq6keytbYeffDyeirwR0HRU85XYlSKsRoAARr9yowvvALgMQMGHEZLH8JGqqdrkYpFUI0AIK1c6nd9h3tbB3Ho+BGaDwAK191uhKlVAjRAAjWzmWQkA6peU5X0n654yB7OBTpNQFKqUM0AIK1c6lt/gmXDuCWRGDs9fbfULrE6WqUUiEiqAAQkckisk5EikXk7jaev0tEVovIchF5X0T6t3huuohsCNymd2bx3cbbCLtXh1/7f0ujrgBPAix61ulKlFIh4pgBICJu4FHgAmAYcJWIDGu12xKgwBgzCngF+H3g2AzgfuBUYDxwv4ikd1753aR8Dfib7BDQcJWQZjuDV7ysncFKKSC4M4DxQLExZpMxphGYAUxtuYMxZp4xpjbw40IgN3D/fGCuMabCGFMJzAUmd07p3WjnMrsN5zMAsNcEaGewUiogmADIAba3+Lkk8NiR3AS8fZzHhqbSpRDXA9IHOF1Jx+QW2M7gRc84XYlSKgQEEwBt9XqaNncUuRYoAP7QnmNF5BYRKRKRovLy8iBK6mY7l0HvUeAK8z7z5s7g0iU21JRSUS2YT7QSoOXYx1ygtPVOInIOcA8wxRjT0J5jjTFPGGMKjDEFWVlZwdbePXxe2L0y/Jt/mjV3Bi/WzmClol0wAVAIDBKRASISC0wDZrXcQUTGAI9jP/zLWjw1BzhPRNIDnb/nBR4LH3vWg7c+PC8Aa0tCGgy/FJa/DA0HnK5GKeWgYwaAMcYL3Ib94F4DzDTGrBKRB0RkSmC3PwDJwMsislREZgWOrQB+jQ2RQuCBwGPho/kK4Eg5AwDbDNRYDatec7oSpZSDPMHsZIyZDcxu9dh9Le6fc5RjnwLCdyayncshJtFOqxwp8sZD1hDbGXzKdU5Xo5RySJj3anaDstX2w9LldrqSztPcGbxjEexa4XQ1SimHaAAcS/layB7qdBWdb9SVdppovTJYqailAXA0tRVwYLc9A4g0iRkwbKqdJrqx9tj7K6UijgbA0ZStsdvs1jNfRIix10PDfu0MVipKaQAcTXlzAETgGQBA/9Mhc5A2AykVpTQAjqZsrZ0Cokf4zV4RlObO4JIv7GynSqmoogFwNGVrbPt/OK4BEKyTrwJ3rF4ZrFQU0gA4mvI1kdv80ywpE4ZOgWUvQlOd09UopbqRBsCRHCiH2r2QFYFDQFsbez3UV8HqN5yuRCnVjTQAjqQs0CYe6WcAAPlnQMYJOk20UlFGA+BIytfabaQOAW1JBMZOh20LbMe3UioqaAAcSdkaiE+D5F5OV9I9Rl8Drhg9C1AqimgAHEnZGjsFRCSPAGopqScMvRiWvaCdwUpFCQ2AthgTGAEUBR3ALRXcYDuDV73udCVKqW6gAdCW6l32gzAaRgC1lH+mnfZ60dNOV6KU6gYaAG2J9CkgjkQExt4A2z+H3aucrkYp1cU0ANoS6ZPAHc3oq+000UV6FqBUpNMAaEv5WkjsaTtGo81h00TXOF2NUqoLaQC0Zc8G6DnY6SqcU3CDnSZ65atOV6KU6kIaAG3ZswF6RtAawO3V7zTbAV74T6crUUp1IQ2A1uoqoXaPnSc/WonAuJtg51K7brBSKiJpALS2p9hue0ZxAIBdMzgmCQqfcroSpVQX0QBobW8gADKjuAkIIL4HjLoCVr5i10ZWSkUcDYDW9m4AlwfS852uxHnjbgJvvV0rQCkVcTQAWtuzwX74u2OcrsR5vUdC3qlQ9JSdHkMpFVE0AFrbWxzdHcCtFdxk/5ts+tDpSpRSnUwDoCW/D/ZuhMwTnK4kdAybComZUPik05UopTqZBkBLVSXga9ARQC3FxMMp02HdbNi33elqlFKdKKgAEJHJIrJORIpF5O42nj9LRBaLiFdELm/1nE9ElgZuszqr8C6xd4PdahPQ4QputNsivTBMqUhyzAAQETfwKHABMAy4SkRaz5K2DbgeeKGNl6gzxowO3KZ0sN6updcAtC0tD066EBY9C031TlejlOokwZwBjAeKjTGbjDGNwAxgassdjDFbjDHLAX8X1Nh99m6AuFRIynK6ktAz/haoq4BVrzldiVKqkwQTADlAy8bfksBjwYoXkSIRWSgil7S1g4jcEtinqLy8vB0v3cn2bLAdwNGyDGR7DDgLsobA54/rkFClIkQwAdDWp2F7PgH6GWMKgKuBv4jIl4bYGGOeMMYUGGMKsrIc/Pa9d6M2/xyJCIy/2c4PVFLkdDVKqU4QTACUAHktfs4FSoN9A2NMaWC7CfgQGNOO+rpPYw3sL9EO4KMZNQ3iesDnjzldiVKqEwQTAIXAIBEZICKxwDQgqNE8IpIuInGB+z2BicDq4y22S+3daLfRPA30scQlw5hvwerXYX/Q3wGUUiHqmAFgjPECtwFzgDXATGPMKhF5QESmAIjIOBEpAb4JPC4izQvKDgWKRGQZMA94yBgTogHQPARUA+CoTr0FjF8vDFMqAniC2ckYMxuY3eqx+1rcL8Q2DbU+bj4wsoM1do/mM4AMvQr4qNLz7ZDQoqfhzB9DbKLTFSmljpNeCdxszwZIzdMPtGBM+J4dErpiptOVKKU6QAOgWcVGyBjodBXhof9EO1Powsd0SKhSYUwDoFnFJp0ELlgiMOFWKF+js4QqFcY0AMCueFVXqWcA7THiG5CUDQsedboSpdRx0gAAqNxstxoAwfPE2ekhiudC2Rqnq1FKHQcNAIAKDYDjMu4m8CTA/EecrkQpdRw0AMC2/4OuA9xeiRkw5lpY/hJU73K6GqVUO2kAgA2AHjkQk+B0JeHntFvB77WTxCmlwooGANgA0Oaf45MxEIZebBeLaah2uhqlVDtoAEAgAAY4XUX4Ov12qK+Cxc85XYlSqh00AOr3Q025TgHREXnjIG8CLPwb+JqcrkYpFSQNAB0C2jnOvAuqtsOKl52uRCkVJA2A5hFAGgAdM+g86DUCPv0z+MN7ZVClooUGwMEA0D6ADhGBM34Ie9bD2jedrkYpFQQNgIpNkNwbYpOcriT8DbsE0gfAp3/SSeKUCgMaABWbtfmns7g9MPEOKF2ik8QpFQY0APQagM41+mp7RvXJH52uRCl1DNEdAI01UL1T2/87kycOJt4OWz6BrfOdrkYpdRTRHQCVW+xWzwA619gb7FTRHz7kdCVKqaOI7gDQIaBdIzbRngVs/gi2LXS6GqXUEWgAgDYBdYWCGyEpS88ClAphGgCJPSE+1elKIk9skp0jaNM82P6F09UopdqgAaDf/rvOuJtswH74/5yuRCnVhigPAL0GoEvFJtnrAjZ+AFs+c7oapVQr0RsA3gaoKrFXrqquM+7b9rqAD36tVwcrFWKiNwD2bQOMNgF1tdhEOPsnsG0BFL/ndDVKqRaiNwCaF4LXM4CuN+Y6SOsP7z+gM4UqFUKCCgARmSwi60SkWETubuP5s0RksYh4ReTyVs9NF5ENgdv0ziq8ww6uA6AB0OU8sfDVn8Ou5bDmDaerUUoFHDMARMQNPApcAAwDrhKRYa122wZcD7zQ6tgM4H7gVGA8cL+IpHe87E5QsQlikuxYddX1Rn4TsobAB7/RVcOUChHBnAGMB4qNMZuMMY3ADGBqyx2MMVuMMcuB1uf35wNzjTEVxphKYC4wuRPq7rjmEUAiTlcSHVxumHQf7C2Gxc86XY1SiuACIAfY3uLnksBjwejIsV2rcjNk5DtdRXQ56ULod7q9Orih2ulqlIp6wQRAW1+Rgx3PF9SxInKLiBSJSFF5eXmQL90Bfj9UbtUO4O4mAuf9BmrK4bOHna5GqagXTACUAHktfs4FSoN8/aCONcY8YYwpMMYUZGV1Q5t8dSn4GrQD2Am5Y2HEN2D+I7A/2F8jpVRXCCYACoFBIjJARGKBacCsIF9/DnCeiKQHOn/PCzzmLB0C6qxJ94HxwQe/dboSpaLaMQPAGOMFbsN+cK8BZhpjVonIAyIyBUBExolICfBN4HERWRU4tgL4NTZECoEHAo85S4eAOis9H079Dix9HkqXOl2NUlHLE8xOxpjZwOxWj93X4n4htnmnrWOfAp7qQI2dr2IzuDzQo82SVXc46yewbAa8/T9w4zs6GkspB0TnlcCVmyGtn13EXDkjPtU2BW1fCCtecboapaJSdAZAxWZt/w8Fo6+FvmNg7i+g4YDT1SgVdaIzACo3a/t/KHC54ILfQ/VO+OSPTlejVNSJvgCorYD6Kj0DCBV54+Hkq2DBI7Cn2OlqlIoq0RcAFToCKOSc8yvwJMBbP9Q1A5TqRtEXAJV6DUDISekF59wHmz+G5TOdrkapqBF9AXDwIrB8R8tQrYy9EXIKYM7PbTOdUqrLRV8AVG6GlD52pSoVOlwuuPhhqKuE9+53uhqlokL0BYAOAQ1dvUfAabfC4n/Blk+drkapiBd9AaBDQEPbV35mA/qN26CxxulqlIpo0RUAjTV2zHnGQKcrUUcSmwRTH7FB/f6vna5GqYgWXQFwcAioBkBIyz8Dxt0Mnz8GW+c7XY1SESvKAmCT3WoAhL5zfglpefDG96Gx1ulqlIpIURYAG+1WAyD0xSXDlEdsaM/9hdPVKBWRoiwANkFSFsT3cLoSFYyBZ8Npt0Hhk7DuHaerUSriRFkAbNZv/+Fm0n3Qa4RtCjpQ5nQ1SkWU6AqAvRsh4wSnq1Dt4YmDbzwJjQdsCOhcQUp1mugJgMZauxi8ngGEn+yhcO4DsOFdOzJIKdUpoicAKrfYrV4EFp7G3wInXQjv/gJKipyuRqmIED0B0DwCKFObgMKSCFzyN+jRB16+XieMU6oTRFEABK4B0HmAwldCOnzzGajeBa9/D/x+pytSKqxFTwDs3QiJmZCQ5nQlqiNyxsL5D8L6d+BTXUZSqY6IngCo2KQjgCLF+Jth5BXwwW9h3dtOV6NU2IqiANBrACKGCEz5K/Q5GV69GcrWOl2RUmEpOgKgqQ72l2gARJKYBJj2gt3OuMouJKOUapfoCIDmIaA6AiiypObAlc/Bvu3w0rfA2+B0RUqFlegIgIOzgOoIoIjTbwJMfRS2fAKzfqBXCivVDh6nC+gWe3UW0Ih28pVQtQ0++A2k9YOv3et0RUqFhaDOAERksoisE5FiEbm7jefjROSlwPOfi0h+4PF8EakTkaWBmzPX8VdsgoQMO45cRaYzfwynXAcf/wGKnnK6GqXCwjHPAETEDTwKnAuUAIUiMssYs7rFbjcBlcaYE0VkGvA74MrAcxuNMaM7ue72qdio3/4jnQhc9Cd7kdibd0FcDxh5udNVKRXSgjkDGA8UG2M2GWMagRnA1Fb7TAWeDdx/BZgkItJ5ZXaQDgGNDu4YuOJf0H8i/Oc7uoaAUscQTADkANtb/FwSeKzNfYwxXqAKyAw8N0BElojIRyJyZltvICK3iEiRiBSVl5e36x9wTI01ULUdMk/s3NdVoSkmAa56EXqPgpnXwaaPnK5IqZAVTAC09U2+9VCLI+2zE+hnjBkD3AW8ICJfWo7LGPOEMabAGFOQlZUVREntsGeD3Wad1Lmvq0JXfA+49lU77PeFK2HTh05XpFRICiYASoC8Fj/nAqVH2kdEPEAqUGGMaTDG7AUwxiwCNgKDO1p0u5Svs9usId36tsphiRlw3Szb9PfClVD8vtMVKRVyggmAQmCQiAwQkVhgGjCr1T6zgOmB+5cDHxhjjIhkBTqREZGBwCBgU+eUHqQ968Dl0T6AaJScBdP/C5mD4MWrYP0cpytSKqQcMwACbfq3AXOANcBMY8wqEXlARKYEdvsnkCkixdimnuahomcBy0VkGbZz+LvGmO6dyL18nf3w98R269uqEJGUCdNnQfYQmHE1LJ/pdEVKhYygLgQzxswGZrd67L4W9+uBb7Zx3KvAqx2ssWPK19olBVX0SsyA6W/aAHjtZqgph9O+73RVSjkusqeC8DbYIaDa/q/ie8A1r8DQKTDn53ZpSV1QRkW5yA6AvRvB+KCnjgBSQEy8XVFs3Ldh/l9h5reg4YDTVSnlmMgOgD3NI4A0AFSAyw0X/i9M/h2smw1PT4aqHU5XpZQjIjsAytcBAj0HdevbVtU1dev7qXYSgQnfhatnQsUWeOJs2PyJ01Up1e0iPADWQnp/e3VoN3nms82c/Kt3mfTHD/nDnLWs3FF1zGP8fp3C2BGDzoVvvwfxafCvqfDZwzqdtIoqkT0ddPn6Lu0AbvL5iXEfytD/LCnhl/9dzWkDM3G54LGPNvHovI2cO6wXD0wdTp9UG0Rl1fW8tngHK0qq2FBWzeY9NfROjef0gT05/cRMzh6cRVqiDlvtFtlD4OYPYNZtMPc+2P4FTPk/O3JIqQgnJsS+8RQUFJiioqKOv5DPCw/2gVO/C+f9uuOv10J1fRP/8+py3lm5i0lDe3HD6fnUNfm45blFjM/P4OkbxhEf46aippEZhdv46/sbcItw61dPZMPuat5asZMmn6F/ZiKDslM4ISuJLXtrWLBxL/vrvcS4hbMHZ3PpmByyUuJYv7ua4rIDjMpN5bJTcjv136ICjIGFf4O590NSFlz6GAw82+mqlAqaiCwyxhS065iIDYC9G+H/ToGpf4Mx13T89QLW7arme/9exNaKWi4e1YePN+yhoqYRgJE5qbxw86mkxMccdsz2ilrueX0lH68vJyXOw+UFuVx3Wj4DeiYdtp/Pb1heso+3lu9k1rJSyqoPLXEY63bR6PPz8wuHcMtZurRllyldCq9+G/YWw+m3wVfvtaOHlApxGgAtrX3LXvjz7Q8gd2zHXw+Ys2oXd85YSnK8h0euGsOpAzOpb/Lx5vKdFG6u4KeTTyIzOa7NY40xrN65n/6ZSSTHHbvlzec3fLG5gnqvj8G9UshOiePOl5by1vKd3HvRUG46YwBFWyt5+rPNHGjwce9FQxncK6VT/p1Rr7EG5twDi56GnoPtkpN5452uSqmj0gBo6ZM/wvsPwN3b7UVAHTRrWSk/fGkpI3NSeeJbY8nu0f3fCr0+P7fPWMLsFbs4ISuJjeU1pCbE4BKoafBx13mDufnMgbhdobMUQ1grfh/+ewdUlcCE78FXfw5xGrIqNB1PAERuJ3D5euiR0ykf/q8uKuEnryyjID+Dp64fF9Q3+K7gcbt4eNoYYtzLWLuzmt9cMoLLTsmhttHHPf9ZwUNvr+X1JTv4yknZjB+Qzpi8dNKTtDP5uJ04CW5dYPsFFv4NVr0O5/8Whl9qh5IqFeYi9wzg8bPtGsDXvd7uQ99cXspzC7YCtm+wcGsFpw3M5MnpBSTGhmZmGmOYtayUZ+ZvYUVJFd7A0NKUOA990xLI75nIZafkMmlINh53ZI/+7RLbC+Gtu2DXchhwNpz/IPQe4XRVSh2kTUDN/H74f7l2kfALHmrXodsrajn3zx+RnRJPn1TbzHNCdjL3fX0Y8THujtXVTeoafSzdvo+VO6rYsa+OHfvqWF6yj937G+ibGs9V4/tx8cl9yW/VCd3MGMOirZUM75tKQmx4/Ju7hd9nF5z/4DdQX2UHF3z1HujR1+nKlNIAOGj3avj7aXDJ32H01UEfZozh+qcLKdpSwXs/OvvguP1I4PX5eW9NGf9euJVPi/cAMKR3CpNH9ObCkX0OdiBvKj/AL95YyWfFexnbP51nbhj3pVFNUa+2wvYxff64XWti/M0w8Q5I6ul0ZSqKaQA0K3oa3rwTfrDYLgsYpFnLSrn9xSXcf/Ewbpg4oGM1hLCSylrmrNrNOyt3UrS1EmPgxOxkRuelMWtpKXEeF5cX5PLcgq2MzE3l2RvH00ND4MsqNsO8B2HFyxCTCKd+x04zrUGgHKAB0Ow/34Xi9+DHG4LurKuqbWLSnz4kJy2B126dGDUjacr21/POql28tdyGwUUj+3Dv14eSnRLPnFW7uO2FxQzr04MHpo6gX0YiaYkxiHaAHq58HXz0O1j5GnjiYcy1NggyIvdLhAo9GgDN/joGsofBtOeD2r2ippEfzVzKxxv2MOu2iQzvm9qx9w9Tfr/B1Sr43lu9m1ufX0yjz86dnxLnYWjfHozLT2dcfgZen2H5jipW7ahiRE4qd0wa9KXXiBrl62H+w7DsJTsN+ZCLYPx3IP8MHTWkupwGAMCBMvjfQXDuA7Zd9hjeWbmTe19fSVVdE/dcOJTrI7jp53iVVNayunQ/2ypq2VZRy7Lt+1hZuh9fYKSRSyA3PZFtFbVcNKoPf7riZOI8bnx+w5vLS9lf7+Wa8f2iJxj274TPH4PFz0Jdpf0yMvYGGPVNOzJNqS6g1wGAncwLIG/CMXe9742V/GvBVkbk9OC5m05laJ+OXzMQiXLTE8lNTzzssZoGL8u27yPG42J43x4kxnp44uONPDh7LZU1jVxRkMcj84opLrMLrny8vpw/XXFydHQo9+gD5/4slT/OAAAOoUlEQVQKvnI3rHgFCv8Bb/8E3r0Xhk2Bk6+yQ0ndkffnp8JL5J0BvPsL++3r7u1HncPlv8tK+cGLS7j+9HzuuWjoYbN6quP32uISfvrKcrx+w4nZyfzwnMHs3l/Pb2evIT8zkb9cOYYTspNIiHFHV1/CzmWw+DlYMdMOIU3KhhGXwfDLIHccuPT3T3WMNgEB/PM8e/XWt+cecZddVfWc9+ePGJiVzCvfPU0vjOpkhVsqKK9u4PzhvQ92pi/ctJfvP7+YvYGJ82LdLtISY8hIiqVnchx90+KZPKI3Zw7KiuwwbqqHDe/aIFg/B3yNkNIHhnzd9hn0nwgevXpbtZ8GgLfBXgB26nfgvN+0uYvfb5j+9BcUbalk9h1nfmlGTtV1yvbX88HaMvbVNVFZ28i+mib21jRSUdPAxvIaquqayEiK5aKRfZgyui9j+6VHdr9BfZUNgdVv2HmHvHUQmwwnfBVOPNdu0/o5XaUKE9oHULrUfqPKO7XNp40x/OOTTXyyYQ8PXjpSP/y7WXaPeKaNb/sDrdHr5+P15by+dAczi7bz3MKt9E2N56JRfTh7cDYF+elhcyV20OJTYdQV9tZYA5s/toGwfg6s+a/dJ2Og7S/IPwP6n65XHatOFVlnAJ/9Feb+wo7/T84++LAxhnnrynj4vQ0sK6ninKHZ/OO6guhqgw4jBxq8vL9mN/9dVspH68tp8hliPS5G56XRu0c8KfEeUuJj6JFgt6kJMUwYkOHIDK1dwhh7bcGmebBxHmydD43V9rm0/nZq6tzxdprzXiPA0/YU5Cq6aBPQjGtg9yq4Y+nBh+qbfNzwdCELNu0lNz2BH3ztRC47JTey25kjSE2Dly+2VDC/eA+FWyqprG2kut7L/rqmgxPeAXhcwjlDe3HNhH6MzEklPsZNnMfFrv31rN1Zzbrd1QzsmcS5w3qFX/D7fbBrBWz9DLYthJJCqN5pn3PFQK9h0Odk6DXSTlDXa7g9u1BRJboDwBg7/v+ESXDZ44GHDHfNXMbrS3fwy4uHc/Wp/fSDP0IYY2jw+tlf30R5dQOzlpby8qKSg6uzHUlB/3TuuWgoY/ql4/cbKmrt/plJseETDMbYNQp2LIKdS6F0iR1lVFd5aJ+Uvna946whkHmivfUcBMm9dcRRhIruPoB926CmHPodav9/6rMt/GfJDu46dzDTT893rjbV6USE+Bg38TFuslPiGd43lbvOG8y8tWWU7qun3uujvslPz+RYhvTuwaDsZOas2sUf567n0r/Np29qPOUHGmjyBabNjvdwQlYyg7KTGZWbyoicVDKT4li+Yx9Lt+2jtKqO3PRE+mcm0j8jid6pdrbYJCfWhhCBtDx7G36JfcwY2F8Ku1faW/k6KF9r58Xy1h061hMP6fmQPgDS+0Nqnu1oTs2x62ckZWtARJGgzgBEZDLwMOAGnjTGPNTq+TjgX8BYYC9wpTFmS+C5nwE3AT7gdmPMnKO9V4eagA6UgzsGEtKYX7yHbz31BZOGZPPYtWMjezSJClpNg5d/frqZTeUH6J2aQJ/UePzGsKm8ho3lB1i7q/pLZxGxHhc5aQns2FdHo9d/2HMpcR6yesSRlRxHemIsFbWN7Kqqp7KmkYL8dC4Zk8O5w3qREOOmusHLnuoGEmM9ZCbHds/ZqN8P1aV2jeM9G6Byi53ErnIz7Nt+qG+hmSsGUnofuiX3sqGQnAVJWZDY0052l5gJ8WkaFiGkS5qARMQNrAfOBUqAQuAqY8zqFvvcCowyxnxXRKYBlxpjrhSRYcCLwHigL/AeMNgY4zvS+3X0OoDNe2r414ItvFS4nb5pCbz+/YmOreClwo8xhtKqelaUVFFR08jInFRO6p1CrMeF32/Ytb+erXtr2b2/nl3769lVVU95dQPl1Q1U1DaSkRRL7x7xJMd7+HBtGaVV9cR5XIhAfdOh8BCBjMRYUhNjbId2vIfc9ASG9U1leN8exHlcbK+oZXtFHV6/oW9aPDlpCaTEx1Db6KW20f4J9UmNp29awvGNkDIG6vfZs+f9pbZZaf8OqN5l+xiqd8GB3Yc3LR1GICHNTm+RkG4DISHN9j/E9bCr8cU135IhNgliUwLbRDvkNSbRdmKHS/NbCOuqJqDxQLExZlPgTWYAU4HVLfaZCvwycP8V4BGxDapTgRnGmAZgs4gUB15vQXuKDMbOqjp+/toK5q0rJ8YtXDiyDz85/yT98FftIiLkpCWQk/bltSBcLqFvWgJ923iuLX6/4YstFby7ajduF2SnxJOZHEtdk4/y6gbKqhuoqm1if30T++u9zF6xixe/2H5cdafEe3AFPkRjPS4G9kxiUK9k+qQmsG1vLRvKqtlWUUdSnJu0hBhSE2PJSIwhIymOzOQkspJHkJU6lqycOFLiPcR53MTHuBAEn7cBc6DcNrHW7IHaPTYUavcidZXENFaR4KvGU7/PnmE07LfXOPiO3h9ziNggiEkIbONtU1VMgg0HT7zduuMC29jDt64Ye+bvjg1sYw495vLY28H77kOPuTwg7sBj7kP3D27l8MfEFbjvCjznOvINCewT2sEWzKdjDtDyt7IEaD3Q/uA+xhiviFQBmYHHF7Y6Nue4qz2K9MRYSvfVc8ekQVxzar/IGRKowpbLJUwYmMmEgZlB7W+MYce+OlaX7sfrN/TLSCQvPRGPWygNrOx2oMFLUpyHpFgPPr9hZ1UdOyrrDl5hDVDb6GVjeQ1vLC2lut5Lz+RYTsxO5pyh2dQ3+aisbWJfbSObyg9QUdN48GwieJmB24mHPZoS56FHQgw+v8EnfmLdjSRSTxK1JFFHstSRQj2JUk+Mv444fz1x/jripZHEpkYSvY0k1DcSL40k0ITbX4fHVBHjt4/FiZc4mojBSyxeYmjCgxc3/jarDBV+BD8CCOawGxhcga19bHvCEIb+z7xuqy2YAGgrwlq3Gx1pn2CORURuAW4B6Nfv+K58jI9x886dZ4bPSA6lWhGRNifeAxjUK4VBgVXbgmWMobbRd8yO6rpGH3sO2DOS8uoGahq8BzvRjTF4XILbJbhcglvsNsYtxLhdxLhd1DX6DjaHVdd7cbvA7ZLD/haNsWdEPmOoMRDrEZrcLupcgjHg8xu8foPX56fJ56fJb0iIcZMU6yYh1oPX5z9Yk9fnxxd4PYNBjA+P8RIrfhJcPuJcPtz4wNeE+L3g9yJ+L8Z4icFHvMsQ7zbEiA+/z4vf58P4veD3YYwPMX7c+PGIwS2GGPHjcYEbP8bvx/h9+P0+XIBbDC78GOPHH3gOYwIf+fbj3SUc/JjHGIyxz9n97M3uA5KWx9B2/V/umGACoATIa/FzLlB6hH1KRMQDpAIVQR6LMeYJ4AmwfQDBFt+afvgrdYiIBDVKKSHWTV5GInkZXw4eFdmC6cIvBAaJyAARiQWmAbNa7TMLmB64fznwgbG9y7OAaSISJyIDgEHAF51TulJKqY445teDQJv+bcAc7DDQp4wxq0TkAaDIGDML+CfwXKCTtwIbEgT2m4ntMPYC3z/aCCCllFLdJ3KuBFZKqSh2PMNA9SoOpZSKUhoASikVpTQAlFIqSmkAKKVUlNIAUEqpKBVyo4BEpBzY2s7DegJ7uqCcrqZ1dy+tu3tp3d3rJGNMuy4XD7mZ0owxWe09RkSK2jv8KRRo3d1L6+5eWnf3EpF2j5/XJiCllIpSGgBKKRWlIiUAnnC6gOOkdXcvrbt7ad3dq911h1wnsFJKqe4RKWcASiml2imsA0BEJovIOhEpFpG7na4nWCKSJyLzRGSNiKwSkTucrilYIuIWkSUi8qbTtbSHiKSJyCsisjbw3/00p2sKhoj8MPA7slJEXhSRkFzqTkSeEpEyEVnZ4rEMEZkrIhsC23Qna2zLEer+Q+D3ZLmI/EdE0pyssS1t1d3iuR+LiBGRnsd6nbANgMBi9Y8CFwDDgKsCi9CHAy/wI2PMUGAC8P0wqv0OYI3TRRyHh4F3jDFDgJMJg3+DiOQAtwMFxpgR2OnYpzlb1RE9A0xu9djdwPvGmEHA+4GfQ80zfLnuucAIY8woYD3ws+4uKgjP8OW6EZE84FxgWzAvErYBQIvF6o0xjUDzYvUhzxiz0xizOHC/Gvth1CVrJXcmEckFLgKedLqW9hCRHsBZ2HUrMMY0GmP2OVtV0DxAQmClvUTaWFEvFBhjPsauBdLSVODZwP1ngUu6taggtFW3MeZdY4w38ONC7EqGIeUI/70B/gz8lDaW3m1LOAdAW4vVh/yHaGsikg+MAT53tpKg/AX7yxXaq3B/2UCgHHg60Hz1pIgkOV3UsRhjdgD/i/02txOoMsa862xV7dLLGLMT7JceINvheo7HjcDbThcRDBGZAuwwxiwL9phwDoCgFpwPZSKSDLwK3GmM2e90PUcjIl8Hyowxi5yu5Th4gFOAvxtjxgA1hGZzxGECbeZTgQFAXyBJRK51tqroISL3YJtrn3e6lmMRkUTgHuC+9hwXzgEQ1ILzoUpEYrAf/s8bY15zup4gTASmiMgWbHPb10Tk386WFLQSoMQY03yW9Qo2EELdOcBmY0y5MaYJeA043eGa2mO3iPQBCGzLHK4naCIyHfg6cI0Jj7HyJ2C/KCwL/I3mAotFpPfRDgrnAAhmsfqQJCKCbY9eY4z5k9P1BMMY8zNjTK4xJh/73/oDY0xYfBs1xuwCtovISYGHJmHXqQ5124AJIpIY+J2ZRBh0XrcwC5geuD8deMPBWoImIpOB/wGmGGNqna4nGMaYFcaYbGNMfuBvtAQ4JfC7f0RhGwCBTprmxerXADONMaucrSpoE4FvYb9FLw3cLnS6qAj3A+B5EVkOjAYedLieYwqcsbwCLAZWYP9eQ/IqVRF5EVgAnCQiJSJyE/AQcK6IbMCOTHnIyRrbcoS6HwFSgLmBv83HHC2yDUeou/2vEx5nN0oppTpb2J4BKKWU6hgNAKWUilIaAEopFaU0AJRSKkppACilVJTSAFBKqSilAaCUUlFKA0AppaLU/wc+tX9KZ8mQNgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = num[0]/(0.5 * N)\n",
    "x = num[1][0:-1] + np.diff(num[1]) * 0.5\n",
    "plt.plot(x,y)\n",
    "plt.hold\n",
    "plt.plot(x, f(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metropolis-Hasting 算法的有效性分析\n",
    "上面的例子展示了 Metropolis-Hasting 算法的有效性。它本质上可以看作是拒绝-接受法的一种扩展。但是对它为何会恰好正确，可能大家还是会有一些疑惑。我们继续用一种直观的方式解释这个现象。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们只考虑有限状态空间，在很多物理实际中，有限已经足够。比如我们考虑一个指定有限区域空间的气体密度变化，我们将这个空间打成细密均匀的网格，每一个网格中的分子数目，本质上就是气体的密度，同时可以看作一种概率状态分布。这是一个有限状态空间。我们用 $m_t(x)$ 表示 $t$ 时刻，空间 $x$ 位置（可以认为是一个空间坐标，或者一个网格编号）的分子数。那么到 $t + 1$ 时刻，$x$ 格子和一个和它相邻的 $y$ 格子之间，会有分子交换，导致 $x$ 格子中气体分子数的变化（在这一过程中，我们只考虑由 $x$ 和 $y$ 两个格子之间交换分子引起的数量变化），这一变化记为：\n",
    "$$\n",
    "\\delta_{xy} = m_t(x) P\\{x \\to y\\} - m_t(y) P\\{y \\to x\\} = m_t(y)P\\{x \\to y\\}\\left(\\frac{m_t(x)}{m_t(y)} - \\frac{P\\{y \\to x\\}}{P\\{x \\to y\\}}\\right).\n",
    "$$\n",
    "这里只用到了质量守恒定律。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据热力学第二定律，在无外界干扰的情况下，经过充分长的时间，气体状态会达到平衡状态，也就是从概率意义上，$\\delta_{xy} = 0$（收支总为零，分布趋于稳定），也即 $\\forall x \\in \\Omega$，有\n",
    "$$\n",
    "\\delta_{xy} = 0 \\Rightarrow \\lim_{t \\to \\infty} \\frac{m_t(x)}{m_t(y)} = \\frac{P\\{y \\to x\\}}{P\\{x \\to y\\}},\n",
    "$$\n",
    "记\n",
    "$$\n",
    "m_e(x) = \\lim_{t \\to \\infty} m_t(x),\n",
    "$$\n",
    "则 $m_e$ 即气体在平衡态下的分布。由 Metropolis-Hasting 算法，\n",
    "$$\n",
    "P\\{x \\to y\\} = g(y | x) h(x, y),\n",
    "$$\n",
    "再由 $g$ 的对称性，有\n",
    "$$\n",
    "\\frac{m_e(x)}{m_e(y)} = \\lim_{t \\to \\infty} \\frac{m_t(x)}{m_t(y)} = \\frac{P\\{y \\to x\\}}{P\\{x \\to y\\}} = \\frac{h(y, x)}{h(x, y)},\n",
    "$$\n",
    "这里最后一项约去了 $g(y | x) = g(x | y)$。而对 $f(x) > f(y)$，$h(y, x) = 1$，此时\n",
    "$$\n",
    "\\frac{m_e(x)}{m_e(y)} = \\frac{1}{h(x, y)} = \\frac{f(x)}{f(y)},\n",
    "$$\n",
    "而对 $f(x) \\leq f(y)$，$h(x, y) = 1$，有\n",
    "$$\n",
    "\\frac{m_e(x)}{m_e(y)} = h(y, x) = \\frac{f(x)}{f(y)},\n",
    "$$\n",
    "由 $x$ 和 $y$ 的任意性，必有 $m_e(x)$ 和 $f(x)$ 成正比例。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "事实上，对平衡态的另一种微观等价描述是：\n",
    "$$\n",
    "f(x)P\\{x \\to y\\} = f(y)P\\{y \\to x\\}, \\quad \\forall x, y \\in \\Omega.\n",
    "$$\n",
    "而 Metropolis-Hasting 算法满足该微观等价性：\n",
    "$$\n",
    "f(x)h(x, y) = f(y)h(y, x).\n",
    "$$\n",
    "这个其实是算法的设计初衷。最后，Metropolis-Hasting 算法产生的马尔可夫链满足：\n",
    "1. 不可约：由 $g$ 的定义知，状态转移可以遍历状态空间 $\\Omega$；\n",
    "2. 非周期：$f$ 至少存在一个最大值 $y^*$，满足其停留的概率是正概率，也即存在可能状态永远停留在 $y^*$。\n",
    "\n",
    "所以该算法产生的马尔可夫链，必存在平稳分布。以上的分析过程并不完全数学上严格，但给出了明确的直观思路。有助于我们理解和掌握此算法，以后辅以概率论技术，也不难严格化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理**  若 $X_t$ 是一个定义在有限状态空间 $\\Omega$ 上的不可约非周期马尔可夫链，其平稳分布为 $\\pi$，设 $\\xi : \\Omega \\to \\mathbb{R}$ 为任意映射，则\n",
    "$$\n",
    "E_\\pi[\\xi] = \\lim_{N \\to \\infty} \\frac{1}{N} \\sum_{i = 1}^N \\xi(x_i)\n",
    "$$\n",
    "由这个定理，我们可以利用 MCMC 方法计算各种随机模型的期望。这里 $\\xi$ 是一个随机变量，$E_\\pi$ 是基于 $\\pi$ 分布的期望。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 统计力学\n",
    "\n",
    "大家可能已经注意到，计算机模拟领域的开拓者们，大部分是物理学家。实际上，计算机模拟这个研究方向的诞生，有三个基本条件：计算机的出现；概率论的成熟；以及物理学（核平）研究的需要。\n",
    "\n",
    "十九世纪以后，人类对物质本质的了解，已经到了原子内部。但即便如此，我们仍然无法直接从分子运动去预测宏观气体的现象，如密度、压力和温度的变化等等。但之前在 MCMC 的理论来源时，我们已经知道，所谓气体的宏观密度，从微观角度看，就是气体分子的空间分布。由于分子数目巨大，我们只要知道其概率分布，就可以对宏观气体密度作出足够好的预测和估算。在类似的思想下，Boltzmann 等人开创了统计力学，建立了一系列描述微观粒子概率分布的方程，比如著名的 Boltzmann 方程：\n",
    "$$\n",
    "\\frac{\\partial}{\\partial t} f(x, v, t) + v \\cdot \\nabla_x f(x, v, t) = Q(f, f) (x, v, t).\n",
    "$$\n",
    "理论上说，对这些方程的统计模拟，可以完全代替宏观模型下的流体方程，如 Euler、Navier-Stikes 等流体力学方程，对流体的宏观性质作出预测和估计。同样的思想，还可以用在电、磁甚至量子力学领域。那么如何正确理解统计物理学模型，如何结合概率论直观给出正确且高效的抽样，并在现代计算机上实现，就称为一项现代科研领域不可缺少的手段。\n",
    "\n",
    "我们这里考虑一个最简单的统计物理模型——**Ising** 模型。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1920 年，德国人 Ernst Ising 在他的博士论文中第一次用统计力学模型解释了铁磁体磁性相变的物理现象。也就是铁磁性物质，存在一个被称为 “居里温度” 的相变点，当温度超过居里温度后，磁性迅速消失，降温至居里温度之下，磁性又会恢复。Ising（实际上已经被我们大大简化）用一个空间网格来模拟铁磁性物质，每个网格代表一个微观铁原子的磁极状态，这个状态这里我们只取两种：北极向上记作 $1$，南极向上记作 $-1$。若大部分铁原子具有相同的磁极状态，那么物质整体就呈现宏观磁性；若原子的磁极是均匀混乱分布的，那么微观磁极互相抵消，宏观物质整体就没有磁性。这里我们假设原子总数是 $N$，那么不同原子的不同磁性选择，会形成 \n",
    "$$\n",
    "K = 2^N\n",
    "$$ 种可能的状态，这是一个巨大的数字，但仍然是有限的。我们这里给出假设：**每一种微观组合状态，都是等概出现的。** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了描述一个系统的混乱程度，Boltzmann 引入了熵（entrpy）的概念：\n",
    "$$\n",
    "S = k_B \\ln K,\n",
    "$$\n",
    "这里 $k_B$ 称为 Boltmann 常数。由热力学第二定律我们知道，熵总是增加的，因此一个平衡系统，熵是极大的（不能再乱了）。我们现在考虑有两个独立的系统，各自的状态数分别是 $K_1$ 和 $K_2$，现在观念上把两个系统合并成一个系统考虑，但实际物理上不考虑它们之间由分子交换引起的新状态，则总系统的状态数为：\n",
    "$$\n",
    "K = K_1 \\cdot K_2,\n",
    "$$\n",
    "而对应地，\n",
    "$$\n",
    "S = k_B \\ln(K_1 \\cdot K_2) = k_B \\ln K_1 + k_B \\ln K_2 = S_1 + S_2,\n",
    "$$\n",
    "即熵有可加性。再由能量守恒定律：\n",
    "$$\n",
    "E = E_1 + E_2 \\Rightarrow E_2 = E - E_1.\n",
    "$$\n",
    "现在假设总能量 $E$ 不随时间变化（但 $E_1$ 和 $E_2$ 会有彼增此减的变化），并将 $S$ 看作 $E_1$ 的函数（也就是看 $E_1$ 如何变化，使得 $S$ 达到最大），则有\n",
    "$$\n",
    "\\left.\n",
    "\\begin{array}{r}\n",
    "\\displaystyle 0 = \\frac{d S}{d E_1} = \\frac{d S_1}{d E_1} + \\frac{d S_2}{d E_2} \\frac{d E_2}{d E_1} \\\\\\\\\n",
    "\\displaystyle \\frac{d E_2}{d E_1} = \\frac{d E}{d E_1} - \\frac{d E_1}{d E_1} = -1\n",
    "\\end{array}\n",
    "\\right\\} \\Rightarrow \\frac{d S_1}{d E_1} - \\frac{d S_2}{d E_2} = 0 \\Rightarrow \\frac{d S_1}{d E_1} = \\frac{d S_2}{d E_2} \\left( = \\frac{d S}{d E}\\right)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过上面的推导，我们的到一个结论：均衡态的 $\\frac{d S}{d E}$ 在系统处处是常数，我们将它记作温度 $T$：\n",
    "$$\n",
    "\\frac{1}{T} := \\frac{d S}{d E}.\n",
    "$$\n",
    "热力学第二定律的另一种表达就是能量总是从高温系统流向低温系统，直至恒温。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去要考虑，对于一个处于平衡态的系统，其中微粒的实际状态服从什么分布。这里首先需要一个基本假设：**只要时间足够长，平衡系统中会经历全部可能的微观状态。**现在我们考虑由 $N$ 个互相不交换粒子的独立系统构成的**正则系综（canonical ensemble）**。它们的总能量是各个系统之和\n",
    "$$\n",
    "E = \\sum_{j = 1}^N E_j,\n",
    "$$\n",
    "总微观状态数是各系统微观状态数之积\n",
    "$$\n",
    "K = \\prod_{j = 1}^N K_j.\n",
    "$$\n",
    "我们令各系统都是保持恒温 $T$ 的，并且系统可能的微观状态共有 $m$ 个，分别是 $\\{\\psi_1, \\psi_2, \\cdots, \\psi_m\\}$，则我们希望能够得到系统处于 $\\psi_i$ 上的概率，也就是概率质量分布。令 $n_i$ 表示处于 $\\psi_i$ 状态的系统数，$E_i$ 表示 $\\psi_i$ 系统所具有的能量，则有\n",
    "$$\n",
    "\\sum_{i = 1}^m n_i = N, \\quad \\sum_{i = 1}^m n_i E_i = E.\n",
    "$$\n",
    "而系综的总微观状态数为\n",
    "$$\n",
    "K = \\frac{N!}{n_1! n_2! \\cdots n_m!}.\n",
    "$$\n",
    "这里分子是系统总数的全排列，分母表示在这个排列下，处于同一状态下的系统应该视作一种情况，因为他们之间交换粒子不产生新的状态。(系综是一种思考模型，考虑足够多的满足恒温 $T$ 条件的系统，它们整体上应该满足什么分布，才符合守恒定律。)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在由最大熵原理，求 $\\{n_i\\}$ 的组合形式，使 $S$ 最大，即\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{rcl}\n",
    "\\displaystyle \\max_{\\{n_i\\}}S & = & k_B \\ln K\\\\\n",
    "\\mathrm{s. t.} && \\sum_{i = 1}^m n_i = N, \\\\\n",
    "&& \\sum_{i = 1}^m n_iE_i = E.\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "由 Lagrange 乘子法，\n",
    "$$\n",
    "\\mathcal{L}\\left(\\{n_i\\}, \\alpha, \\beta\\right) = k_B \\ln K - \\alpha \\sum_{i = 1}^m n_i - \\beta \\sum_{i = 1}^m n_i E_i,\n",
    "$$\n",
    "其中 $\\alpha$ 和 $\\beta$ 是对应的 Lagrange 乘子。于是\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial n_i} & = & k_B\\frac{\\partial}{\\partial n_i} \\ln K - \\alpha \\frac{\\partial}{\\partial n_i} \\sum_{j = 1}^m n_j - \\beta \\frac{\\partial}{\\partial n_i}\\sum_{j = 1}^m n_j E_i \\\\\n",
    "& = & k_B \\frac{\\partial}{\\partial n_i} \\ln K - \\alpha - \\beta E_i = 0.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "注意这里第一行的第二项当 $i = j$ 时是 $1$，其余都是 $0$；第三项当 $i = j$ 时是 $E_i$，其余是 $0$。而对第二行第一项，由 Stirling 公式，\n",
    "$$\n",
    "\\ln n! \\approx n \\ln n - n,\n",
    "$$\n",
    "故\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\ln \\frac{N!}{n_1!\\cdots n_m!} & = & \\ln N! - \\left(\\sum_{i = 1}^m \\ln n_i\\right) \\\\\n",
    "& = & N \\ln N - N - \\sum_{i} n_i \\ln n_i + \\sum_{i} n_i\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "注意到 \n",
    "$$\n",
    "N = \\sum_i n_i,\n",
    "$$\n",
    "于是\n",
    "$$\n",
    "\\ln \\frac{N!}{n_1!\\cdots n_m!} = N \\ln N - \\sum_i n_i \\ln n_i.\n",
    "$$\n",
    "进而\n",
    "$$\n",
    "\\frac{\\partial \\ln K}{\\partial n_i} = \\frac{\\partial}{\\partial n_i}\\left(N\\ln N - \\sum_j n_j \\ln n_j\\right),\n",
    "$$\n",
    "这里第一项是常数，求导后为零；第二项根据求导乘法展开，只有 $i = j$ 才非零，故\n",
    "$$\n",
    "\\frac{\\partial \\ln K}{\\partial n_i} = -\\frac{\\partial}{\\partial n_i} n_i \\ln n_i = -\\ln n_i - 1,\n",
    "$$\n",
    "综合有，\n",
    "$$\n",
    "-k_B \\ln n_i - k_B - \\alpha - \\beta E_i = 0,\n",
    "$$\n",
    "所以\n",
    "$$\n",
    "n_i = e^{-1 - \\frac{\\alpha}{k_B} - \\beta \\frac{E_i}{k_B}},\n",
    "$$\n",
    "和\n",
    "$$\n",
    "e^{-\\frac{\\beta}{k_B}E_i}\n",
    "$$\n",
    "成正比。因此，系统在各状态 $\\Omega = \\{\\psi_1, \\psi_2, \\cdots, \\psi_N\\}$ 上的概率密度函数（这里其实还是概率质量）为\n",
    "$$\n",
    "p_i = \\frac{e^{-\\frac{\\beta}{k_B}}E_i}{Z},\n",
    "$$\n",
    "其中 $Z$ 是归一化因子，满足\n",
    "$$\n",
    "Z = \\sum_{i = 1}^m e^{-\\frac{\\beta}{k_B}}E_i,\n",
    "$$\n",
    "使得\n",
    "$$\n",
    "\\sum_{i \\in \\Omega}p_i = 1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下去我们估计 $\\beta$，为此我们首先修改一下模型中熵的定义：\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{rcl}\n",
    "\\displaystyle \\max_{n_i} S & = & -k_B \\sum_{i = 1}^m p_i \\log p_i \\\\\n",
    "\\mathrm{s. t.} && \\sum_{i = 1}^m p_i = 1\\\\\n",
    "&& \\sum_{i = 1}^m p_i E_i = \\bar{E}\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "这里采用的熵定义是信息论大佬 Shannon 引入的信息熵：\n",
    "$$\n",
    "S = - \\sum_{i = 1}^m p_i \\ln p_i,\n",
    "$$\n",
    "$p_i$ 是某概率分布。它和 Boltzmann 引入的熵是等价的（差一个常数因子）。我们仍然使用 Lagrange 乘子法来分析：\n",
    "$$\n",
    "\\mathcal{L}\\left(\\{p_i\\}, \\alpha, \\beta\\right) = -k_b \\sum_{i = 1}^m p_i \\ln p_i - \\alpha \\sum_{i = 1}^m p_i - \\beta \\sum_{i = 1}^m p_i E_i,\n",
    "$$\n",
    "则\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial p_i} = \\sum_{i = 1}^m (-k_B \\ln p_i - k_B - \\alpha - \\beta E_i) = 0,\n",
    "$$\n",
    "对任意的 $\\{p_i\\}$ 均成立，因此\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "&& -k_B \\ln p_i - k_B - \\alpha - \\beta E_i  =  0, \\\\\n",
    "&\\Rightarrow& -k_B \\ln p_i  =  k_B + \\alpha + \\beta E_i, \\\\\n",
    "&\\Rightarrow& -k_B p_i \\ln p_i  =  k_B p_i + \\alpha p_i + \\beta E_i p_i, \\\\\n",
    "&\\Rightarrow& -k_B \\sum_{i = 1}^m p_i \\ln p_i  =  \\sum_{i = 1}^m (k_B p_i + \\alpha p_i + \\beta E_i p_i), \\\\\n",
    "&& =  k_B + \\alpha + \\beta \\bar{E}, \\\\\n",
    "&\\Rightarrow& S = k_B + \\alpha + \\beta \\bar{E}.\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "由温度定义，在平衡态，由\n",
    "$$\n",
    "\\beta = \\frac{\\partial S}{\\partial \\bar{E}} = \\frac{1}{T}.\n",
    "$$\n",
    "代入 $p_i$ 的表达式，有\n",
    "$$\n",
    "p_i = \\frac{1}{Z}e^{-\\frac{E_i}{k_BT}}, \\quad Z = \\sum_{i = 1}^m e^{-\\frac{E_i}{k_BT}}.\n",
    "$$\n",
    "这个分布，被称为 Boltzmann 分布。它为我们微观模拟提供了抽样依据。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们终于可以回到 Ising 模型。铁磁体的总磁矩量是一个宏观量，但它事实上是由单个原子自旋产生的微观磁场的叠加。现在假设目标铁磁体处于热平衡状态，而各原子的微观量服从 Boltzmann 分布，而宏观量本质上就是这些微观量的一种加权叠加，或称期望：\n",
    "\n",
    "$$\n",
    "<A> = \\sum_{x \\in \\Omega} A_x\\frac{e^{-\\frac{E_x}{k_BT}}}{Z},\n",
    "$$\n",
    "\n",
    "其中 $A_x$ 代表在 $x$ 位置的微观量，$<A>$ 称为 Boltzmann 期望，对应我们观念中的宏观量。在有了 Metropolis 算法之后，我们可以不用纠结如何求归一化系数 $Z$，于是一切都可以实际抽样和模拟。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实际算例\n",
    "\n",
    "我们来考虑一个简化版的 Ising 模型。用一个 $N \\times N$ 的等距二维网格来模拟一块铁磁体。每个网格交叉点代表一个铁原子，它具有 $\\pm 1$ 两种自选状态之一。即\n",
    "$$\n",
    "s(i, j) = \\pm 1, \\quad \\Omega = \\{S | S = s(i, j),\\quad i, j = 1, 2, \\cdots, N\\},\n",
    "$$\n",
    "而总磁矩为\n",
    "$$\n",
    "M(S) = \\frac{1}{N^2}\\sum_{j = 1}^N \\sum_{i = 1}^N s(i, j).\n",
    "$$\n",
    "显然，若 $s(i, j)$ 是均匀分布的，那么 $M(S) \\approx 0$.\n",
    "\n",
    "对于一个实际模型，其边界，也就是铁磁体如何和外界交换信息是重要的。在本问题中，我们试图回避这个问题，或者，我们的理想情况就是一块无限边界的铁磁体，这样我们可以对有限区域的铁磁体采用**周期边界条件**：\n",
    "$$\n",
    "\\begin{array}{rclrcl}\n",
    "s(0, j) &=& s(N, j), & s(N + 1, j) &=& s(1, j),\\\\\n",
    "s(i, 0) &=& s(i, N), & s(i, N + 1) &=& s(i, 1),\\\\\n",
    "s(0, 0) &=& s(N, N), & s(0, N + 1) &=& s(N, 1),\\\\\n",
    "s(N + 1, 0) &=& s(1, N), & s(N + 1, N + 1) &=& s(1, 1).\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "这样相当于模拟一块无限但是结构循环出现的物质。我们也可以假设不存在外界，或者物质与外界彻底隔绝，表现为我们关心的量，比如这里的 $M$，延边界法向的变化率为零，也即\n",
    "$$\n",
    "\\frac{\\partial M}{\\partial \\vec{n}} = 0.\n",
    "$$\n",
    "这种边界也称为**自由边界条件**。在实际抽样时，直接当外部格点不存在就行了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在 $\\forall (i, j)$，定义邻域为\n",
    "$$\n",
    "\\Lambda (i, j) = \\{(i - 1, j), (i, j - 1), (i, j + 1), (i + 1, j)\\},\n",
    "$$\n",
    "也就是和 $(i, j)$ 直接相邻的格点。而能量则由经验公式\n",
    "$$\n",
    "E(S) = -\\frac{1}{2} J \\sum_{j = 1}^N \\sum_{i = 1}^N \\sum_{\\mu \\in \\Lambda (i, j)} s(\\mu) s(i, j),\n",
    "$$\n",
    "决定，这里 $J > 0$ 是常数。而由之前分析，当系统处于热平衡态时，状态为 $s \\in \\Omega$ 的概率服从 Boltzmann 分布，即\n",
    "$$\n",
    "P(S) = \\frac{1}{Z}e^{-\\frac{E(S)}{k_BT}}, \\quad S \\in \\Omega.\n",
    "$$\n",
    "因此我们不难用随机模拟来估计总磁矩量\n",
    "\n",
    "$$\n",
    "<M> = \\frac{1}{Z} \\sum_{S \\in \\Omega} M(S) e^{-\\beta E(S)}, \\quad \\beta = \\frac{1}{k_BT}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们用 Metropolis-Hasting 算法来模拟抽样，这里可以采用最小邻域（变化）思想，即每一次的建议改变，都是在现有的一个状态 $S$ 中随机抽取一个格点，然后对其微磁矩进行反转，得到一个新的状态 $S'$，然后根据\n",
    "$$\n",
    "h(S, S') = \\min \\left\\{1, \\frac{P(S')}{P(S)}\\right\\} = \\min\\left\\{1, \\frac{e^{-\\beta E(S')}/Z}{e^{-\\beta E(SS)}/Z}\\right\\} = \\min \\left\\{1, e^{-\\beta[E(S')-E(S)]}\\right\\}\n",
    "$$\n",
    "来判定是否接受建议改变。下面是模拟的代码。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1 -1 -1]\n",
      " [ 1 -1  1]\n",
      " [-1  1 -1]]\n"
     ]
    }
   ],
   "source": [
    "N = 3\n",
    "S = np.random.randint(0, 2, (N, N)) * 2 - 1\n",
    "print(S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Energy(S):\n",
    "    E = S[0][0] * S[0][1] + S[0][0] * S[1][0]\n",
    "    E = E + S[N - 1][0] * S[N - 2][0] + S[N - 1][0] * S[N - 1][1]\n",
    "    E = E + S[0][N - 1] * S[0][N - 2] + S[0][N - 1] * S[1][N - 1]\n",
    "    E = E + S[N - 1][N - 1] * S[N - 1][N - 2] + S[N - 1][N - 1] * S[N - 2][N - 1]\n",
    "    for i in range(1, N - 1):\n",
    "        E = E + S[i][0] * S[i - 1][0] + S[i][0] * S[i + 1][0] + S[i][0] * S[i][1]\n",
    "        E = E + S[i][N - 1] * S[i - 1][N - 1] + S[i][N - 1] * S[i + 1][N - 1] + S[i][N - 1] * S[i][N - 2]\n",
    "        for j in range(1, N - 1):\n",
    "            E = E + S[0][j] * S[0][j - 1] + S[0][j] * S[0][j + 1] + S[0][j] * S[1][j]\n",
    "            E = E + S[N - 1][j] * S[N - 1][j - 1] + S[N - 1][j] * S[N - 1][j + 1] + S[N - 1][j] * S[N - 2][j]\n",
    "            E = E + S[i][j] * S[i - 1][j] + S[i][j] * S[i + 1][j] + S[i][j] * S[i][j + 1] + S[i][j] * S[i][j - 1]\n",
    "    E = -E * 0.5\n",
    "    return E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Sample(warmup, trials, S, T):\n",
    "    M = np.zeros(trials)\n",
    "    beta = 1 / T\n",
    "    Es = Energy(S)\n",
    "    for t in range(warmup + trials):\n",
    "        i = np.random.randint(N)\n",
    "        j = np.random.randint(N)\n",
    "        Y = S.copy()\n",
    "        Y[i][j] = -Y[i][j]\n",
    "        Ey = Energy(Y)\n",
    "        H = np.min([1, np.exp(-beta * (Ey - Es))])\n",
    "        if np.random.rand() < H:\n",
    "            S = Y.copy()\n",
    "            Es = Ey\n",
    "        if t >= warmup:\n",
    "            M[t - warmup] = np.sum(S)\n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = Sample(100000, 100000, S, 100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  206.,     0.,  1941.,     0.,  7340.,     0., 16328.,     0.,\n",
       "        24029.,     0., 24266.,     0., 16530.,     0.,  7245.,     0.,\n",
       "         1897.,   218.]),\n",
       " array([-9., -8., -7., -6., -5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,\n",
       "         4.,  5.,  6.,  7.,  8.,  9.]),\n",
       " <a list of 18 Patch objects>)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAEh5JREFUeJzt3X+M3PV95/Hn6+wQtU16OMEQDpwzrayqTqUSYhFX6CSu9MCQqian5GT+ACvlzmkO7hqpleKkUolCK5E7JVXRJVSksWKqNAQl4bASp46PixSdVCgLpfwIybGlbtjaBaemhBOnRk7f/WM+m079mfWOd5ed2fJ8SKP5zvv7+X73Pd8Z72u/P2acqkKSpGH/YtINSJKmj+EgSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkzvpJN7BU55xzTm3evHnSbUjSmvLwww9/t6o2LjZu0XBIsgm4C3gT8PfAnVX1u0k+DPwn4Hgb+qGqOtiW+SBwI/AD4L9W1aFW3wH8LrAO+P2quq3VLwLuBt4APAJcX1XfP11fmzdvZmZmZrH2JUlDkvzlOOPGOax0Evi1qvppYDtwU5Ktbd7vVNXF7TYfDFuBXcBbgB3AJ5OsS7IO+ARwNbAVuG5oPR9t69oCvMAgWCRJE7JoOFTVsap6pE2/BDwFXHCaRXYCd1fV31XVXwCzwKXtNltVz7S9gruBnUkC/Dzwhbb8fuDapT4hSdLyndEJ6SSbgbcCD7bSzUkeS7IvyYZWuwB4dmixuVZbqP5G4G+r6uQp9VE/f0+SmSQzx48fHzVEkrQCxg6HJK8Dvgi8v6q+B9wB/CRwMXAM+Nj80BGL1xLqfbHqzqraVlXbNm5c9HyKJGmJxrpaKclrGATDZ6vqSwBV9dzQ/E8BX24P54BNQ4tfCBxt06Pq3wXOTrK+7T0Mj5ckTcCiew7tnMCngaeq6uND9fOHhr0TeKJNHwB2JXltuwppC/AnwEPAliQXJTmLwUnrAzX434a+DryrLb8buG95T0uStBzj7DlcBlwPPJ7k0Vb7EIOrjS5mcAjoCPBegKp6Msk9wDcZXOl0U1X9ACDJzcAhBpey7quqJ9v6PgDcneS3gD9lEEaSpAnJWv1vQrdt21Z+zkGSzkySh6tq22Lj/PoMSVJnzX59hvTP2ea9X1n2Oo7c9o4V6ESvVu45SJI67jlIQ/yLXRpwz0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1DEcJEkdw0GS1Fk0HJJsSvL1JE8leTLJr7b6G5IcTvJ0u9/Q6klye5LZJI8luWRoXbvb+KeT7B6qvy3J422Z25PklXiykqTxjLPncBL4tar6aWA7cFOSrcBe4P6q2gLc3x4DXA1sabc9wB0wCBPgFuDtwKXALfOB0sbsGVpux/KfmiRpqRYNh6o6VlWPtOmXgKeAC4CdwP42bD9wbZveCdxVAw8AZyc5H7gKOFxVJ6rqBeAwsKPN+/Gq+uOqKuCuoXVJkibgjM45JNkMvBV4EDivqo7BIECAc9uwC4Bnhxaba7XT1edG1CVJEzJ2OCR5HfBF4P1V9b3TDR1RqyXUR/WwJ8lMkpnjx48v1rIkaYnGCockr2EQDJ+tqi+18nPtkBDt/vlWnwM2DS1+IXB0kfqFI+qdqrqzqrZV1baNGzeO07okaQnGuVopwKeBp6rq40OzDgDzVxztBu4bqt/QrlraDrzYDjsdAq5MsqGdiL4SONTmvZRke/tZNwytS5I0AevHGHMZcD3weJJHW+1DwG3APUluBL4DvLvNOwhcA8wCLwPvAaiqE0luBR5q4z5SVSfa9PuAzwA/Any13SRJE7JoOFTV/2H0eQGAK0aML+CmBda1D9g3oj4D/MxivUiSVoefkJYkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkddZPugFJ02vz3q8sex1HbnvHCnSi1WY4aCr4S0iaLh5WkiR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1Fg2HJPuSPJ/kiaHah5P8VZJH2+2aoXkfTDKb5NtJrhqq72i12SR7h+oXJXkwydNJPp/krJV8gpKkMzfOnsNngB0j6r9TVRe320GAJFuBXcBb2jKfTLIuyTrgE8DVwFbgujYW4KNtXVuAF4Abl/OEJEnLt2g4VNU3gBNjrm8ncHdV/V1V/QUwC1zabrNV9UxVfR+4G9iZJMDPA19oy+8Hrj3D5yBJWmHLOedwc5LH2mGnDa12AfDs0Ji5Vluo/kbgb6vq5Cl1SdIELTUc7gB+ErgYOAZ8rNUzYmwtoT5Skj1JZpLMHD9+/Mw6liSNbUnhUFXPVdUPqurvgU8xOGwEg7/8Nw0NvRA4epr6d4Gzk6w/pb7Qz72zqrZV1baNGzcupXVJ0hiWFA5Jzh96+E5g/kqmA8CuJK9NchGwBfgT4CFgS7sy6SwGJ60PVFUBXwfe1ZbfDdy3lJ4kSStn0a/sTvI54HLgnCRzwC3A5UkuZnAI6AjwXoCqejLJPcA3gZPATVX1g7aem4FDwDpgX1U92X7EB4C7k/wW8KfAp1fs2UmSlmTRcKiq60aUF/wFXlW/Dfz2iPpB4OCI+jP842EpSdIU8BPSkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6iwaDkn2JXk+yRNDtTckOZzk6Xa/odWT5PYks0keS3LJ0DK72/ink+weqr8tyeNtmduTZKWfpCTpzIyz5/AZYMcptb3A/VW1Bbi/PQa4GtjSbnuAO2AQJsAtwNuBS4Fb5gOljdkztNypP0uStMoWDYeq+gZw4pTyTmB/m94PXDtUv6sGHgDOTnI+cBVwuKpOVNULwGFgR5v341X1x1VVwF1D65IkTchSzzmcV1XHANr9ua1+AfDs0Li5VjtdfW5EfaQke5LMJJk5fvz4EluXJC1mpU9IjzpfUEuoj1RVd1bVtqratnHjxiW2KElazFLD4bl2SIh2/3yrzwGbhsZdCBxdpH7hiLokaYKWGg4HgPkrjnYD9w3Vb2hXLW0HXmyHnQ4BVybZ0E5EXwkcavNeSrK9XaV0w9C6JEkTsn6xAUk+B1wOnJNkjsFVR7cB9yS5EfgO8O42/CBwDTALvAy8B6CqTiS5FXiojftIVc2f5H4fgyuifgT4artJkiZo0XCoqusWmHXFiLEF3LTAevYB+0bUZ4CfWawPSdLq8RPSkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqTOov8TnP7527z3K8ta/sht71ihTiRNC8NB0tRb7h8w4B8xZ8rDSpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEgSeoYDpKkjuEgSeosKxySHEnyeJJHk8y02huSHE7ydLvf0OpJcnuS2SSPJblkaD272/ink+xe3lOSJC3XSuw5/NuquriqtrXHe4H7q2oLcH97DHA1sKXd9gB3wCBMgFuAtwOXArfMB4okaTJeicNKO4H9bXo/cO1Q/a4aeAA4O8n5wFXA4ao6UVUvAIeBHa9AX5KkMS03HAr4WpKHk+xptfOq6hhAuz+31S8Anh1adq7VFqpLkiZkuf9N6GVVdTTJucDhJN86zdiMqNVp6v0KBgG0B+DNb37zmfYqSRrTsvYcqupou38euJfBOYPn2uEi2v3zbfgcsGlo8QuBo6epj/p5d1bVtqratnHjxuW0Lkk6jSWHQ5IfS/L6+WngSuAJ4AAwf8XRbuC+Nn0AuKFdtbQdeLEddjoEXJlkQzsRfWWrSZImZDmHlc4D7k0yv54/rKo/SvIQcE+SG4HvAO9u4w8C1wCzwMvAewCq6kSSW4GH2riPVNWJZfQlSVqmJYdDVT0D/OyI+t8AV4yoF3DTAuvaB+xbai+SpJXlJ6QlSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSR3DQZLUMRwkSZ31k27g1Wzz3q8sex1HbnvHCnQiSf+U4SBJY3i1/THnYSVJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUsdwkCR1DAdJUmdqwiHJjiTfTjKbZO+k+5GkV7OpCIck64BPAFcDW4HrkmydbFeS9Oo1Ld/KeikwW1XPACS5G9gJfPOV+GGvtm9XlDQd1tLvnqnYcwAuAJ4dejzXapKkCUhVTboHkrwbuKqq/mN7fD1waVX9l1PG7QH2tIc/BXx7VRs9M+cA3510E2NaK73a58paK33C2ul1LfT5r6tq42KDpuWw0hywaejxhcDRUwdV1Z3AnavV1HIkmamqbZPuYxxrpVf7XFlrpU9YO72ulT7HMS2HlR4CtiS5KMlZwC7gwIR7kqRXranYc6iqk0luBg4B64B9VfXkhNuSpFetqQgHgKo6CBycdB8raE0c/mrWSq/2ubLWSp+wdnpdK30uaipOSEuSpsu0nHOQJE0Rw2GFJPl8kkfb7UiSRxcYdyTJ423czGr32Xr4cJK/Gur3mgXGTfQrTZL89yTfSvJYknuTnL3AuIls08W2T5LXtvfFbJIHk2xerd6GetiU5OtJnkryZJJfHTHm8iQvDr0ffnO1+xzq5bSvZQZub9v0sSSXTKDHnxraVo8m+V6S958yZmq26ZJVlbcVvgEfA35zgXlHgHMm3N+HgV9fZMw64M+BnwDOAv4M2LrKfV4JrG/THwU+Oi3bdJztA/xn4Pfa9C7g8xN4rc8HLmnTrwf+74g+Lwe+vNq9LeW1BK4BvgoE2A48OOF+1wF/zeCzA1O5TZd6c89hhSUJ8B+Az026l2X64VeaVNX3gfmvNFk1VfW1qjrZHj7A4PMv02Kc7bMT2N+mvwBc0d4fq6aqjlXVI236JeAp1va3D+wE7qqBB4Czk5w/wX6uAP68qv5ygj28IgyHlfdvgOeq6ukF5hfwtSQPt098T8rNbbd8X5INI+ZP21ea/DKDvxhHmcQ2HWf7/HBMC7kXgTeuSncjtMNabwUeHDH755L8WZKvJnnLqjb2Ty32Wk7b+3IXC/8hOC3bdEmm5lLWtSDJ/wLeNGLWb1TVfW36Ok6/13BZVR1Nci5wOMm3quobq9krcAdwK4N/iLcyOAz2y6euYsSyK35p2zjbNMlvACeBzy6wmlXZpqcYZ/usyjYcR5LXAV8E3l9V3ztl9iMMDov8v3b+6X8CW1a7x2ax13KatulZwC8BHxwxe5q26ZIYDmegqn7hdPOTrAf+PfC206zjaLt/Psm9DA5PrPgvssV6nZfkU8CXR8wa6ytNlmuMbbob+EXgimoHc0esY1W26SnG2T7zY+bae+NfAide4b46SV7DIBg+W1VfOnX+cFhU1cEkn0xyTlWt+ncEjfFarsr7ckxXA49U1XOnzpimbbpUHlZaWb8AfKuq5kbNTPJjSV4/P83ghOsTq9jffB/Dx2jfuUAPE/9KkyQ7gA8Av1RVLy8wZlLbdJztcwDY3abfBfzvhQLuldLOcXwaeKqqPr7AmDfNnwtJcimD3wt/s3pd/rCPcV7LA8AN7aql7cCLVXVslVudt+BRgmnZpsvhnsPK6o4/JvlXwO9X1TXAecC97T2zHvjDqvqjVe8S/luSixnsjh8B3ntqrzUdX2nyP4DXMji8APBAVf3KNGzThbZPko8AM1V1gMEv5T9IMstgj2HXK93XCJcB1wOP5x8vr/4Q8GaAqvo9BsH1viQngf8P7FrtEGtGvpZJfmWo14MMrliaBV4G3jOBPknyo8C/o/3babXhPqdlmy6Zn5CWJHU8rCRJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqTOPwCGFZwGKYP7QQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(M, bins = 18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里有一个坑，Python 默认的数组复制是传址的，比如："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1. 1.]\n",
      "[1. 9. 1. 1.]\n",
      "[1. 9. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "a = np.ones(4)\n",
    "print(a)\n",
    "b = a\n",
    "b[1] = 9\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意到这里 a 和 b 一起被修改。正确的复制操作是对 list："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 3, 4]\n",
      "[1, 2, 9, 4]\n"
     ]
    }
   ],
   "source": [
    "a = [1, 2, 3, 4]\n",
    "b = a[:]\n",
    "b[2] = 9\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对 array："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4]\n",
      "[1 2 9 4]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([1, 2, 3, 4])\n",
    "b = a.copy()\n",
    "b[2] = 9\n",
    "print(a)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ising(T):    # T is the tempure.  +-+-+   Ising area.\n",
    "    nTrials = 100000    # for stat.          |3|4|   N = 2.\n",
    "    startup = 100000    # for steady.        +-+-+   free boundary.\n",
    "    d = np.zeros(nTrials)   # results.       |1|2| \n",
    "    beta = 1/T             #                +-+-+ \n",
    "    M = 0\n",
    "    s = 2 * np.random.randint(0, 2, 5) - 1   # randomly set 1 or -1.\n",
    "    Es = s[1] * s[2] + s[3] * s[4]        \n",
    "    Es = Es + s[1] * s[3] + s[2] * s[4]   \n",
    "    Es = -Es                              \n",
    "    for t in range(startup + nTrials):        \n",
    "        k = np.random.randint(1, 5)  # randomly pick 1 ~ 4.\n",
    "        #print(\"k, S = \", k, s, Es)\n",
    "        y = s        # y is the suggestion dist in s' neighbour. \n",
    "        y[k] = -s[k]  # randomly flip the status on one point.\n",
    "        #print(\"s = \",  s)\n",
    "        Ey = y[1] * y[2] + y[3] * y[4]\n",
    "        Ey = Ey + y[1] * y[3] + y[2] * y[4]\n",
    "        Ey = -Ey\n",
    "        #print(\"k, Y = \", k, y, Ey)\n",
    "        h = np.min([1, np.exp(-beta * (Ey - Es))])\n",
    "        if (np.random.rand() < h):\n",
    "            s = y\n",
    "            Es = Ey\n",
    "            #print(s, Es)\n",
    "        if (t >= startup):\n",
    "            Ms = s[1] + s[2] + s[3] + s[4]  \n",
    "            M = M + Ms                      \n",
    "            d[t - startup] = Ms\n",
    "        \n",
    "    \n",
    "#    x = -4:1:4;\n",
    "    plt.hist(d);\n",
    "#    M = (M / nTrials) / 4;\n",
    "#    D = sum(d.^2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFKpJREFUeJzt3W2sXeWZ3vH/NQYSNJnUEA4pxU6NWmsaQieGuMYVX1JIwZBRTKoggdpgZag8jaBKpLQNJFKZhCARjSZUqAkVU1xMmwlBeRFWMHVcQhRFCi+HxLw4DuWU0HACxU4NhAiVCHL3w36sbPnZ9nmzvQ/4/5OW9l73etZa9wbOuc5e69mbVBWSJA37g3E3IElafAwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdY4ZdwPzddJJJ9WKFSvG3YYkvaE8/PDDv6qqiZnGvWHDYcWKFUxOTo67DUl6Q0nyv2czzstKkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqSO4SBJ6hgOkqTOG/YT0tJiteLqu8d27qdv+ODYzq03F985SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6M4ZDkrcmeTDJI0l2Jvlcq9+W5OdJdrRlVasnyU1JppI8muSsoWNtSPJkWzYM1d+X5LG2z01JcjherCRpdmbz9RmvAudW1W+SHAv8MMk9bdu/rapv7Df+QmBlW84GbgbOTnIicC2wGijg4SRbquqFNmYjcD+wFVgH3IMkaSxmfOdQA79pq8e2pQ6yy3rg9rbf/cDSJKcAFwDbq2pvC4TtwLq27e1V9aOqKuB24OIFvCZJ0gLN6p5DkiVJdgC7GfyCf6Btur5dOroxyVta7VTgmaHdp1vtYPXpEfVRfWxMMplkcs+ePbNpXZI0D7MKh6p6vapWAcuANUnOAK4B/gHwj4ATgU+34aPuF9Q86qP6uKWqVlfV6omJidm0LkmahznNVqqqF4HvA+uq6rl26ehV4L8Aa9qwaWD50G7LgGdnqC8bUZckjclsZitNJFnanh8PfAD4WbtXQJtZdDHweNtlC3B5m7W0Fnipqp4DtgHnJzkhyQnA+cC2tu3lJGvbsS4H7jq0L1OSNBezma10CrA5yRIGYXJnVX0nyfeSTDC4LLQD+Fdt/FbgImAKeAX4GEBV7U1yHfBQG/f5qtrbnn8cuA04nsEsJWcqSdIYzRgOVfUocOaI+rkHGF/AlQfYtgnYNKI+CZwxUy+SpCPDT0hLkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpYzhIkjqGgySpM2M4JHlrkgeTPJJkZ5LPtfppSR5I8mSSryc5rtXf0tan2vYVQ8e6ptWfSHLBUH1dq00lufrQv0xJ0lzM5p3Dq8C5VfVeYBWwLsla4IvAjVW1EngBuKKNvwJ4oar+PnBjG0eS04FLgfcA64CvJFmSZAnwZeBC4HTgsjZWkjQmM4ZDDfymrR7blgLOBb7R6puBi9vz9W2dtv28JGn1O6rq1ar6OTAFrGnLVFU9VVW/Be5oYyVJYzKrew7tL/wdwG5gO/C/gBer6rU2ZBo4tT0/FXgGoG1/CXjHcH2/fQ5UlySNyazCoaper6pVwDIGf+m/e9Sw9pgDbJtrvZNkY5LJJJN79uyZuXFJ0rzMabZSVb0IfB9YCyxNckzbtAx4tj2fBpYDtO1/C9g7XN9vnwPVR53/lqpaXVWrJyYm5tK6JGkOZjNbaSLJ0vb8eOADwC7gPuAjbdgG4K72fEtbp23/XlVVq1/aZjOdBqwEHgQeAla22U/HMbhpveVQvDhJ0vwcM/MQTgE2t1lFfwDcWVXfSfJT4I4kXwB+Atzaxt8K/NckUwzeMVwKUFU7k9wJ/BR4Dbiyql4HSHIVsA1YAmyqqp2H7BVKkuZsxnCoqkeBM0fUn2Jw/2H/+v8DLjnAsa4Hrh9R3wpsnUW/kqQjwE9IS5I6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6hoMkqWM4SJI6M4ZDkuVJ7kuyK8nOJJ9o9b9I8sskO9py0dA+1ySZSvJEkguG6utabSrJ1UP105I8kOTJJF9PctyhfqGSpNmbzTuH14BPVdW7gbXAlUlOb9turKpVbdkK0LZdCrwHWAd8JcmSJEuALwMXAqcDlw0d54vtWCuBF4ArDtHrkyTNw4zhUFXPVdWP2/OXgV3AqQfZZT1wR1W9WlU/B6aANW2Zqqqnquq3wB3A+iQBzgW+0fbfDFw83xckSVq4Y+YyOMkK4EzgAeAc4KoklwOTDN5dvMAgOO4f2m2a34fJM/vVzwbeAbxYVa+NGC/pDWDF1XeP5bxP3/DBsZz3aDDrcEjyNuCbwCer6tdJbgauA6o9/hXwZ0BG7F6MfpdSBxk/qoeNwEaAd73rXbNtXfjDK2luZjVbKcmxDILhq1X1LYCqer6qXq+q3wF/zeCyEQz+8l8+tPsy4NmD1H8FLE1yzH71TlXdUlWrq2r1xMTEbFqXJM3DbGYrBbgV2FVVXxqqnzI07MPA4+35FuDSJG9JchqwEngQeAhY2WYmHcfgpvWWqirgPuAjbf8NwF0Le1mSpIWYzWWlc4CPAo8l2dFqn2Ew22gVg0tATwN/DlBVO5PcCfyUwUynK6vqdYAkVwHbgCXApqra2Y73aeCOJF8AfsIgjCRJYzJjOFTVDxl9X2DrQfa5Hrh+RH3rqP2q6il+f1lKkjRmfkJaktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktSZMRySLE9yX5JdSXYm+USrn5hke5In2+MJrZ4kNyWZSvJokrOGjrWhjX8yyYah+vuSPNb2uSnJqP9ntSTpCJnNO4fXgE9V1buBtcCVSU4HrgburaqVwL1tHeBCYGVbNgI3wyBMgGuBs4E1wLX7AqWN2Ti037qFvzRJ0nzNGA5V9VxV/bg9fxnYBZwKrAc2t2GbgYvb8/XA7TVwP7A0ySnABcD2qtpbVS8A24F1bdvbq+pHVVXA7UPHkiSNwZzuOSRZAZwJPAC8s6qeg0GAACe3YacCzwztNt1qB6tPj6iPOv/GJJNJJvfs2TOX1iVJczDrcEjyNuCbwCer6tcHGzqiVvOo98WqW6pqdVWtnpiYmKllSdI8zSockhzLIBi+WlXfauXn2yUh2uPuVp8Glg/tvgx4dob6shF1SdKYzGa2UoBbgV1V9aWhTVuAfTOONgB3DdUvb7OW1gIvtctO24Dzk5zQbkSfD2xr215Osrad6/KhY0mSxuCYWYw5B/go8FiSHa32GeAG4M4kVwC/AC5p27YCFwFTwCvAxwCqam+S64CH2rjPV9Xe9vzjwG3A8cA9bZEkjcmM4VBVP2T0fQGA80aML+DKAxxrE7BpRH0SOGOmXiRJR4afkJYkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVJnxnBIsinJ7iSPD9X+Iskvk+xoy0VD265JMpXkiSQXDNXXtdpUkquH6qcleSDJk0m+nuS4Q/kCJUlzN5t3DrcB60bUb6yqVW3ZCpDkdOBS4D1tn68kWZJkCfBl4ELgdOCyNhbgi+1YK4EXgCsW8oIkSQs3YzhU1Q+AvbM83nrgjqp6tap+DkwBa9oyVVVPVdVvgTuA9UkCnAt8o+2/Gbh4jq9BknSILeSew1VJHm2XnU5otVOBZ4bGTLfagervAF6sqtf2q0uSxmi+4XAz8PeAVcBzwF+1ekaMrXnUR0qyMclkksk9e/bMrWNJ0qzNKxyq6vmqer2qfgf8NYPLRjD4y3/50NBlwLMHqf8KWJrkmP3qBzrvLVW1uqpWT0xMzKd1SdIszCsckpwytPphYN9Mpi3ApUnekuQ0YCXwIPAQsLLNTDqOwU3rLVVVwH3AR9r+G4C75tOTJOnQOWamAUm+BrwfOCnJNHAt8P4kqxhcAnoa+HOAqtqZ5E7gp8BrwJVV9Xo7zlXANmAJsKmqdrZTfBq4I8kXgJ8Atx6yVydJmpcZw6GqLhtRPuAv8Kq6Hrh+RH0rsHVE/Sl+f1lKkrQI+AlpSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVJnxnBIsinJ7iSPD9VOTLI9yZPt8YRWT5KbkkwleTTJWUP7bGjjn0yyYaj+viSPtX1uSpJD/SIlSXMzm3cOtwHr9qtdDdxbVSuBe9s6wIXAyrZsBG6GQZgA1wJnA2uAa/cFShuzcWi//c8lSTrCZgyHqvoBsHe/8npgc3u+Gbh4qH57DdwPLE1yCnABsL2q9lbVC8B2YF3b9vaq+lFVFXD70LEkSWMy33sO76yq5wDa48mtfirwzNC46VY7WH16RF2SNEaH+ob0qPsFNY/66IMnG5NMJpncs2fPPFuUJM1kvuHwfLskRHvc3erTwPKhccuAZ2eoLxtRH6mqbqmq1VW1emJiYp6tS5JmMt9w2ALsm3G0AbhrqH55m7W0FnipXXbaBpyf5IR2I/p8YFvb9nKStW2W0uVDx5IkjckxMw1I8jXg/cBJSaYZzDq6AbgzyRXAL4BL2vCtwEXAFPAK8DGAqtqb5DrgoTbu81W17yb3xxnMiDoeuKctkqQxmjEcquqyA2w6b8TYAq48wHE2AZtG1CeBM2bqQ5J05PgJaUlSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSx3CQJHUMB0lSZ8b/n8Ob0Yqr7x7LeZ++4YNjOa+kQ+/N/nvEdw6SpI7hIEnqGA6SpM6CwiHJ00keS7IjyWSrnZhke5In2+MJrZ4kNyWZSvJokrOGjrOhjX8yyYaFvSRJ0kIdincO/6SqVlXV6rZ+NXBvVa0E7m3rABcCK9uyEbgZBmECXAucDawBrt0XKJKk8Tgcl5XWA5vb883AxUP122vgfmBpklOAC4DtVbW3ql4AtgPrDkNfkqRZWmg4FPDdJA8n2dhq76yq5wDa48mtfirwzNC+0612oHonycYkk0km9+zZs8DWJUkHstDPOZxTVc8mORnYnuRnBxmbEbU6SL0vVt0C3AKwevXqkWMkSQu3oHcOVfVse9wNfJvBPYPn2+Ui2uPuNnwaWD60+zLg2YPUJUljMu9wSPKHSf5o33PgfOBxYAuwb8bRBuCu9nwLcHmbtbQWeKlddtoGnJ/khHYj+vxWkySNyUIuK70T+HaSfcf5m6r670keAu5McgXwC+CSNn4rcBEwBbwCfAygqvYmuQ54qI37fFXtXUBfkqQFmnc4VNVTwHtH1P8vcN6IegFXHuBYm4BN8+1FknRo+QlpSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdQwHSVLHcJAkdRZNOCRZl+SJJFNJrh53P5J0NFsU4ZBkCfBl4ELgdOCyJKePtytJOnotinAA1gBTVfVUVf0WuANYP+aeJOmotVjC4VTgmaH16VaTJI1BqmrcPZDkEuCCqvqXbf2jwJqq+tf7jdsIbGyrfww8Mc9TngT8ap77Hk72NTf2NTf2NTdv1r7+blVNzDTomAWc4FCaBpYPrS8Dnt1/UFXdAtyy0JMlmayq1Qs9zqFmX3NjX3NjX3NztPe1WC4rPQSsTHJakuOAS4EtY+5Jko5ai+KdQ1W9luQqYBuwBNhUVTvH3JYkHbUWRTgAVNVWYOsROt2CL00dJvY1N/Y1N/Y1N0d1X4vihrQkaXFZLPccJEmLyFEfDkn+TZJKctK4ewFIcl2SR5PsSPLdJH9n3D0BJPnLJD9rvX07ydJx9wSDadBJdib5XZKxzyxZjF8Dk2RTkt1JHh93L8OSLE9yX5Jd7d/hJ8bdE0CStyZ5MMkjra/PjbunYUmWJPlJku8czvMc1eGQZDnwT4FfjLuXIX9ZVX9SVauA7wD/ftwNNduBM6rqT4D/CVwz5n72eRz4Z8APxt3IIv4amNuAdeNuYoTXgE9V1buBtcCVi+Sf16vAuVX1XmAVsC7J2jH3NOwTwK7DfZKjOhyAG4F/ByyaGy9V9euh1T9kkfRWVd+tqtfa6v0MPosydlW1q6rm+2HIQ21Rfg1MVf0A2DvuPvZXVc9V1Y/b85cZ/MIb+zcj1MBv2uqxbVkUP4dJlgEfBP7z4T7XURsOST4E/LKqHhl3L/tLcn2SZ4B/zuJ55zDsz4B7xt3EIuTXwMxTkhXAmcAD4+1koF262QHsBrZX1aLoC/gPDP6g/d3hPtGimcp6OCT5H8DfHrHps8BngPOPbEcDB+urqu6qqs8Cn01yDXAVcO1i6KuN+SyDywFfPRI9zbavRSIjaoviL87FLMnbgG8Cn9zvnfPYVNXrwKp2b+3bSc6oqrHes0nyp8Duqno4yfsP9/ne1OFQVR8YVU/yD4HTgEeSwOASyY+TrKmq/zOuvkb4G+BujlA4zNRXkg3AnwLn1RGcAz2Hf17jNquvgdHvJTmWQTB8taq+Ne5+9ldVLyb5PoN7NuO+oX8O8KEkFwFvBd6e5L9V1b84HCc7Ki8rVdVjVXVyVa2oqhUMfqjPOhLBMJMkK4dWPwT8bFy9DEuyDvg08KGqemXc/SxSfg3MHGTwl9mtwK6q+tK4+9knycS+2XhJjgc+wCL4Oayqa6pqWfuddSnwvcMVDHCUhsMid0OSx5M8yuCy16KY3gf8R+CPgO1tmu1/GndDAEk+nGQa+MfA3Um2jauXdsN+39fA7ALuXAxfA5Pka8CPgD9OMp3kinH31JwDfBQ4t/03taP9VTxupwD3tZ/Bhxjcczis00YXIz8hLUnq+M5BktQxHCRJHcNBktQxHCRJHcNBktQxHCRJHcNBktQxHCRJnf8PlOh+PZgMMzoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Ising(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
